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A  numerical  method  has  been  developed  for  the  unsteady  Navier-Stokes  equa¬ 
tions  of  incompressible  flow  in  three  dimensions.  The  momentum  equations,  com¬ 
bined  with  a  pressure  correction  equation,  are  solved  employing  a  non-st aggered  grid 
which  results  in  a  simpler  formulation  compared  to  the  classical  approach  of  using 
staggered  meshes.  The  momentum  equations  are  solved  explicitly  using  a  finite  vol¬ 
ume  algorithm,  while  the  pressure  Poisson  equation  is  discretized  using  the  Galerkin 
finite  element  method  and  implicitly  solved.  The  grid  is  formed  with  hybrid  (pris¬ 
matic/tetrahedral)  elements.  Equation  adaptation  is  utilized,  the  Navier-Stokes 
equations  are  solved  in  the  prismatic  region,  which  includes  the  viscous  region,  and 
the  Euler  equations  are  solved  in  the  tetrahedral  region,  which  is  inviscid.  Adap¬ 
tive  local  grid  refinement  of  the  prism  and  tetrahedral  cells  is  employed  in  order  to 
optimize  the  mesh  to  the  flow  solution.  Validation  of  the  algorithm  is  performed 
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using  experimental  and  other  numerical  data,  and  demonstrate  the  accuracy  and 
robustness  of  the  unsteady  three-dimensional  method. 

An  additional  study  is  conducted  using  an  existing  two-dimensional  incom¬ 
pressible  Navier-Stokes  solver.  The  two-dimensional  solver  is  applied  to  predict  the 
hydrodynamic  forces  on  a  circular  cylinder  due  to  reversing  flows. 
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Chapter  1 


Introduction 


Incompressible  flows  are  frequently  encountered  in  engineering  applications  in  the 
fields  of  aerodynamics,  hydraulics,  meteorology  and  hydrodynamics.  Applications 
of  incompressible  flows  are  also  found  in  specialized  topics  such  as  stratified  flows, 
turbulence,  biological  fluid  mechanics,  analysis  of  cloud  droplets  and  ice  crystals, 
sediment  transport,  rotating  flows  and  lubrication  theory. 

1.1  Incompressible  versus  Compressible  Flow  Equations 

The  incompressible  flow  equations  are  a  subset  of  the  compressible  flow  equations. 
Incompressible  flows  could  be  solved  using  the  general  compressible  flow  equations. 
For  compressible  flow  solvers,  the  time  step  used  to  march  the  solution  in  time  is 
proportional  to  the  inverse  of  the  speed  of  sound.  As  the  incompressible  flow  limit 
is  approached,  the  theoretical  speed  of  sound  approaches  infinity.  This  results  in 
the  time  step  of  compressible  flow  solvers  approaching  zero.  As  a  consequence,  nu¬ 
merical  simulation  of  incompressible  flows  with  compressible  flow  solvers  becomes 
prohibitively  slow  and  therefore  computationally  expensive.  Implicit  schemes  are 
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often  used  to  increase  the  time  step  limitation  of  explicit  time  marching  methods. 
However,  the  increase  in  time  step  is  practically  less  than  a  factor  of  5-10  [7]  due  to 
the  growth  in  truncation  error.  Therefore,  for  both  explicit  and  implicit  schemes, 
solution  of  the  compressible  governing  equations  is  prohibitively  expensive  for  in¬ 
compressible  flows.  It  is  advantageous  to  employ  the  incompressible  Navier-Stokes 
equations  rather  than  the  compressible  equations.  The  incompressible  Navier-Stokes 
equations  appear  simpler  than  the  compressible  form  of  the  Navier-Stokes  equations 
when  isothermal  flow  is  assumed  because  the  energy  equation  is  decoupled  from  the 
momentum  and  continuity  equations.  If  the  temperature  is  not  of  interest,  the 
energy  equation  does  not  need  to  be  solved  for  incompressible  flows. 

1.2  Issues  Related  to  the  Numerical  Solution  of  the  In¬ 
compressible  Flow  Equations 

A  few  major  issues,  related  to  solving  numerically  the  incompressible  flow  equations, 
must  be  addressed.  The  infinite  theoretical  speed  of  sound  in  incompressible  flows 
requires  special  treatment  in  the  solution  procedure  for  the  pressure  field.  The 
solution  of  the  pressure  field  is  often  decoupled  from  the  solution  of  the  velocity 
field.  Several  distinct  methods  have  been  developed  to  solve  incompressible  flow 
equations  which  address  the  decoupling  of  the  pressure  and  velocity  fields.  A  second 
issue  of  importance  is  the  use  of  staggered  or  non-staggered  grids.  The  location  of 
the  pressure  and  velocity  unknowns  have  significant  consequences  on  the  solver’s 
convergence  and  solution.  A  third  major  issue  is  the  topology  of  the  computational 
grid  employed.  Structured  grids  are  easier  to  use  and  the  flow  solver  is  simpler  to 
develop.  Application  of  structured  grid  solvers  are  limited  as  it  is  very  difficult  to 
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develop  a  structured  grid  about  a  complex  three-  dimensional  surface.  Unstructured 
grids  have  the  capability  to  cover  a  complex  shape,  but  unstructured  flow  solvers 
typically  place  substantial  demands  on  computer  storage  requirements  as  well  as 
being  more  difficult  to  develop.  Historical  development  of  these  issues  is  presented 
as  well  as  discussion  of  the  methodology  chosen  for  the  present  work. 

1.2.1  Pressure/ Velocity  Decoupling 

During  the  past  two  decades  a  significant  number  of  numerical  algorithms  have 
been  developed  for  the  solution  of  the  incompressible  Navier- Stokes  equations.  In 
compressible  flows,  the  conservation  of  mass  is  given  as  a  partial  differential  equation 
for  the  temporal  variation  of  density.  In  fact,  the  set  of  five  equations:  continuity, 
x,  y,  z-momentum  and  energy,  can  be  solved  explicitly  for  the  unknown  density, 

u, v,w  velocities,  and  energy,  since  each  equation  contains  the  temporal  variation 
of  the  unknown.  The  pressure  is  then  determined  from  the  energy  equation  using 
an  equation  of  state.  The  four  incompressible  equations  have  the  unknowns:  u, 

v,  w  velocities  and  pressure.  Lack  of  a  pressure  term  in  the  continuity  equation 
makes  solution  of  the  momentum  equations  with  the  divergence-free  constraint  more 
difficult.  In  the  case  of  incompressible  flows,  conservation  of  mass  acts  as  a  constraint 
condition  that  the  velocity  field  needs  to  satisfy.  Furthermore,  spatial  discretization 
of  the  pressure  and  velocity  may  produce  oscillatory  solutions. 

The  three  major  categories  of  algorithms  used  to  solve  incompressible  flows 
are  shown  in  Figure  1.1.  Historically,  the  first  category  (vorticity  based  formulations) 
calculates  the  values  of  the  stream  function  and  vorticity,  and  determines  the  velocity 
and  pressure  fields  afterwards.  The  other  two  categories,  artificial  compressibility 
and  pressure  correction,  both  use  the  primitive  variables,  pressure  and  velocity,  as 
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Figure  1.1:  Major  categories  of  incompressible  flow  algorithms, 
where  ^  :  stream  function 
u  :  vorticity 
V  :  velocity 
(f)  :  potential 


unknowns  but  are  significantly  different  in  design. 

The  vorticity  based  formulation  completely  decouples  the  velocity  and  pres¬ 
sure  calculations.  In  two-dimensional  calculations,  the  governing  equations  are  for¬ 
mulated  in  terms  of  a  stream  function  and  a  vorticity  [1,  2,  6,  26,  33,  40,  59, 
110,  119,  137].  Direct  extension  of  this  method  to  three  dimensions  is  not  possi¬ 
ble.  However,  different  formulations  have  been  used  in  three-dimensions,  such  as 
the  vorticity- velocity  approach  [96],  the  vector-potential  approach  [132],  and  the 
vector  stream  function  method  [116].  The  selection  of  boundary  conditions  is  often 
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quite  challenging  for  these  methods.  The  methods  also  seem  computationally  ex¬ 
pensive  since  multiple  Laplace  or  Poisson  equations  need  to  be  solved  for  each  time 
iteration. 

The  artificial  compressibility  approach  uses  compressible-like  governing  equa¬ 
tions.  The  artificial  compressibility  approach  in  two-dimensions  [28,  30,  111,  140] 
and  in  three-dimensions  [73,  121,  123,  146]  adds  a  time  derivative  of  the  pressure  to 
the  continuity  equation  and  the  incompressible  flow  field  is  treated  as  compressible 
during  the  transient  stage.  Time  accuracy  of  the  simulation  is  usually  not  preserved, 
but  recent  progress  has  been  made  in  the  development  of  time  accurate  algorithms 
in  two  dimensions  [16,  74,  82,  84,  95,  112,  113,  140]  and  three  dimensions  [85,  114]. 
Time  accurate  algorithms  usually  require  additional  iterations  in  pseudo-time  to 
converge  the  solution  at  each  time  step. 

The  third  class  of  algorithms,  called  pressure  correction  methods,  use  a  Pois¬ 
son  equation  for  the  pressure  field  [50,  99].  The  usual  computational  procedure 
is  to  assume  an  initial  pressure  field,  and  then  an  iterative  process  is  defined  until 
the  continuity  equation  is  satisfied.  The  present  work  uses  a  pressure  correction 
formulation. 

1.2.2  Staggered  versus  Non-Staggered  Grids 

A  major  consideration  of  pressure  correction  methods  is  whether  the  grid  is  stag¬ 
gered  or  non-staggered.  A  non-staggered  grid  has  velocity  components  and  pressure 
collocated  with  each  other.  In  the  finite  element  context  using  similar  approxima¬ 
tion  functions  for  the  velocity  and  pressure  is  analogous  to  the  use  of  non-staggered 
grids.  The  use  of  the  non-staggered  grid  can  produce  oscillations  in  the  pressure 
field.  In  the  finite  element  context  this  is  related  to  the  inf-sup  condition  [24,  25]. 
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In  order  to  reject  these  oscillatory  modes  staggered  grids,  or  different  approxima¬ 
tion  functions  for  finite  element  methods,  have  been  employed  by  several  pressure 
correction  type  algorithms  in  two-dimensions  [8,  15,  20,  21,  23,  27,  32,  39,  42, 
50,  55,  56,  68,  71,  88,  89,  100,  103,  128,  129,  134,  135]  and  in  three  dimensions 
[19,  41,  47,  48,  51,  52,  60,  67,  69,  70,  72,  87,  127].  Use  of  these  staggered  grids  with 
local  embedding  adaptive  grids  complicates  the  treatment  of  grid  interfaces. 

On  the  other  hand,  employment  of  non-staggered  grids  requires  some  special 
treatment  of  the  pressure  correction  formulation  to  couple  the  velocity  and  pres¬ 
sure  and  eliminate  the  oscillations.  One  method  modifies  the  continuity  equation 
in  the  derivation  of  the  pressure  correction  equations  to  include  a  slight  compress¬ 
ibility  effect  in  two-dimensions  [58,  80,  124,  141]  and  three  dimensions  [10,  49]. 
A  second  method  uses  the  momentum  interpolation  method  in  two- dimensions 
[13,  14,  29,  78,  83,  86,  102,  107]  and  in  three-dimensions  [79,  122,  147].  The  mo¬ 
mentum  interpolation  method  assumes  the  velocities  at  any  control  volume  face  are 
driven  by  the  pressure  difference  between  the  two  adjacent  nodes  on  either  side  of 
the  face  in  question,  which  is  similar  to  the  basic  concept  of  staggered  grids.  A 
linear  interpolation  is  used  to  evaluate  cell  face  convective  and  diffusive  terms  from 
the  cell  centered  convective  and  diffusive  terms.  A  third  method  requires  dissipa¬ 
tion  in  the  algorithms  in  two-dimensions  [64,  81,  91,  141]  and  three-dimensions 
[9,  106,  125,  126,  130].  Stability  of  both  approaches  with  high  Reynolds  number 
flows  is  an  important  issue.  Staggered  and  non-staggered  grids  will  be  discussed 
in  greater  detail  in  the  section  on  Numerical  Issues.  The  present  work  uses  a  non- 
staggered  grid  and  artificial  dissipation  to  eliminate  oscillations. 
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1.2.3  Unstructured  versus  Structured  Meshes 


Simulation  of  flows  around  three-dimensional  bodies  is  a  major  issue  in  computa¬ 
tional  fluid  mechanics.  The  generation  of  a  structured  hexahedral  grid  which  covers 
a  complex  three  dimensional  body  has  been  proven  to  be  very  difficult  [11].  Un¬ 
structured  tetrahedral  grids  provide  flexibility  in  three-dimensional  grid  generation 
since  they  can  cover  complicated  topologies  more  easily  compared  to  the  hexahe¬ 
dral  meshes.  Unstructured  grids  have  been  used  extensively  for  three-dimensional 
inviscid  flow  modeling  of  compressible  flows.  Substantial  progress  is  currently  being 
made  in  large  scale  viscous  flow  computations  on  complex  three-dimensional  geome¬ 
tries  using  compressible  flow  solvers.  However,  there  has  been  very  little  analogous 
work  in  two-dimensions  [81,  143]  or  in  three  dimensions  [75,  77]  using  incompressible 
viscous  flow  solvers. 

For  viscous  flow  solvers,  it  is  critical  that  the  grid  have  clustered  points  in 
the  normal  direction  to  any  no- slip  boundaries.  This  usually  results  in  high  aspect 
ratio  grid  cells  in  the  viscous  region.  It  is  difficult  to  generate  high  aspect  ratio  cells 
using  an  unstructured  grid  while  minimizing  the  number  of  cells.  Unstructured  grids 
using  tetrahedral  cells  also  typically  place  substantial  demands  on  computer  storage 
requirements  because  of  indirect  addressing  and  arbitrariness  of  nodal  distribution 
which  requires  the  need  for  extensive  connectivity  information.  It  is  therefore  im¬ 
perative  that  a  numerical  scheme  be  storage  efficient.  Most  existing  unstructured 
Navier- Stokes  methods  need  a  supercomputer  to  compute  viscous  flows  on  complex 
three-dimensional  geometries.  One  of  the  prime  objectives  of  the  present  work  has 
been  to  implement  a  scheme  that  could  fit  viscous  simulations  in  a  workstation. 

The  mesh  methodology  proposed  by  [66,  90]  provides  an  alternative  ap¬ 
proach  to  generate  suitable  grids  for  Navier-Stokes  computations  by  using  prismatic 
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elements.  The  prisms  are  unstructured  in  the  lateral  direction,  while  they  are  struc¬ 
tured  in  the  normal- to- surface  direction.  The  prisms  permit  the  use  of  sufficient 
grid  clustering  in  the  normal  direction  as  well  as  flexibility  in  covering  complex  sur¬ 
faces  using  their  triangular  faces.  Prismatic  data  structures  also  reduce  the  storage 
requirements  considerably  compared  to  tetrahedral  data  structures. 

The  hybrid  mesh  [63,  138]  incorporates  both  prismatic  and  tetrahedral  ele¬ 
ments.  Viscous  effects  are  dominant  in  a  relatively  small  region  of  the  flow  domain, 
in  close  proximity  to  the  no-slip  wall.  Inviscid  flow  features  that  are  dominant 
away  from  the  body  do  not  exhibit  the  directionality  that  is  characteristic  of  vis¬ 
cous  stresses.  If  prismatic  cells  are  used  in  the  viscous  portion  of  the  flow,  they 
will  fill  only  a  small  fraction  of  the  computational  domain.  Since  unstructured  grids 
have  the  capability  to  fill  a  given  volume  with  fewer  cells,  the  inviscid  portion  of 
the  domain  is  tessellated  with  tetrahedra  which  are  generated  starting  from  the 
triangulation  of  the  outermost  surface  of  the  prismatic  region. 

1.3  Present  Research 

The  primary  goal  of  this  dissertation  is  to  develop  and  evaluate  a  three-  dimen¬ 
sional  incompressible  Navier-Stokes  method  using  hybrid  (prismatic/tetrahedral) 
grids.  The  flow  domain  is  discretized  with  prismatic  and  tetrahedral  elements  using 
a  non-staggered  grid.  The  use  of  the  hybrid  mesh  enables  the  solver  to  compute 
accurate  solutions  with  relatively  less  memory  requirements  than  a  fully  unstruc¬ 
tured  mesh.  The  use  of  prismatic  cells  to  discretize  the  full  Navier-Stokes  equations 
and  tetrahedra  cells  to  discretize  the  Euler  equations  renders  the  solver  equation- 
adaptive.  The  present  work  is  the  first  procedure  which  computes  three-dimensional 
incompressible  viscous  flows  with  hybrid  elements. 
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A  secondary  goal  is  to  apply  grid  adaptation  to  the  present  method.  A  hybrid 
grid  adaptation  scheme  implementing  local  mesh  refinement  provides  optimum  grids 
for  viscous  flow  computations.  Grid  refinement  is  a  dual  adaptation  scheme  that 
couples  a  three-dimensional  isotropic  division  of  tetrahedra  and  two-dimensional 
directional  division  of  prisms.  The  advantage  of  using  local  adaptation  is  examined; 
namely,  yielding  accurate  results  with  reduced  computing  resources. 

A  tertiary  goal  is  to  study  the  forces  on  a  circular  cylinder  in  an  oscillating 
freestream  flow.  This  study  uses  an  existing  two-dimensional  incompressible  Navier- 
Stokes  solver. 


1.4  Overview  of  Dissertation 

A  brief  summary  of  the  contents  of  each  chapter  follows. 

Chapter  2  details  the  pressure  correction  numerical  method  used  in  the 
present  work.  The  Navier-Stokes  equations  are  non-dimensionalized  and  discretized 
using  a  non-staggered  grid.  Equations  for  the  pressure  correction  method  are  de¬ 
rived  and  methodology  described.  The  numerical  method  was  developed  using  the 
ideas  and  formulations  of  many  different  authors.  The  synthesis  of  these  different 
techniques  resulted  in  a  numerical  method  which  is  unique. 

Chapter  3  describes  the  adaptive  hybrid  (prismatic/tetrahedral)  grid.  Equa¬ 
tion  adaptation  as  well  as  the  spatial  adaptation  methodology  is  detailed  and  in¬ 
cludes  a  description  of  the  flow  feature  detection  algorithm. 

Chapter  4  details  the  integration  scheme  developed  in  the  present  work.  The 
concept  of  the  dual  cell  used  in  the  edge-based  integration  is  described,  as  well  as  the 
finite  volume  discretization  of  the  momentum  equations.  Finite  element  discretiza¬ 
tion  of  the  pressure  Poisson  equation  is  described  and  includes  the  construction  of  the 
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global  matrix.  Initial  conditions  are  described  and  special  considerations  required 
for  boundary  conditions  on  unstructured  grids  is  presented  as  well  as  the  artificial 
dissipation  methodology,  time  step  calculation,  and  solver  computer  requirements. 

Chapter  5  details  the  major  numerical  issues  encountered  in  the  development 
of  this  three-dimensional  incompressible  adaptive  Navier-Stokes  solver.  Pressure 
correction  formulations  are  discussed,  and  the  importance  of  diagonal  dominance 
and  boundary  conditions  for  the  Poisson  equation  are  emphasized.  The  consequence 
of  non-staggered  versus  staggered  grids  is  presented  as  well  as  an  explanation  of  the 
presence  of  a  non-zero  discrete  divergence.  Accuracy  of  spatial  discretization  is 
determined  by  testing  with  analytic  field  functions. 

Chapter  6  presents  the  validation  test  cases  of  the  adaptive  hybrid  solver. 
Flow  over  a  flat  plate  is  used  to  demonstrate  the  stability  of  the  solver  at  high 
Reynolds  number.  A  quasi-two-dimensional  cylinder,  three-dimensional  cubic  cav¬ 
ity  and  sphere  are  used  to  demonstrate  the  validity  of  the  all-prismatic  solver  for 
steady  flows.  An  impulsive  start  of  flow  around  the  cylinder  is  used  to  validate 
time  accuracy  of  the  solver.  The  cubic  driven  cavity  case  is  used  to  illustrate  the 
usefulness  of  prismatic  adaptation,  and  the  sphere  is  used  to  demonstrate  hybrid 
adaptation.  A  special  case  of  flow  over  tandem  spheres  is  used  to  demonstrate  the 
suitability  of  hybrid  adaptation. 

Chapter  7  presents  a  summary  of  the  present  work,  and  includes  the  contri¬ 
butions  of  this  research,  conclusions,  and  recommendations  for  further  research. 

Appendix  E  presents  the  study  of  the  hydrodynamic  forces  on  a  circular 
cylinder  in  an  oscillating  freestream  flow  using  an  existing  two-dimensional  solver. 
The  simulations  are  conducted  for  two  oscillating  flow  regimes  and  the  lift  and  in-line 
force  coefficients  are  compared  against  both  numerical  and  experimental  results. 
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Chapter  2 


Pressure  Correction  Numerical 

Method 


The  laws  of  conservation  of  mass  and  momentum  for  fluid  flow  are  expressed  in 
integral  form  for  a  three-dimensional  domain.  The  three-dimensional  Navier-Stokes 
equations  are  then  derived  and  expressed  in  non-dimensional  form.  The  semi¬ 
discrete  form  of  the  governing  equations  are  presented  and  the  pressure  correction 
formulation  is  described. 

2.1  3-D  Navier-Stokes  Equations 

Consider  a  volume  f 2  enclosed  by  a  boundary  surface  d£l  with  h  being  the  outward 
normal  on  the  boundary  surface.  The  velocity  of  the  fluid  is  given  by  V,  and  the 
fluid  in  the  control  volume  is  assumed  to  be  non-reacting  and  isothermal. 

The  law  of  conservation  of  mass  requires  the  time  rate  of  decrease  of  mass  of 
the  fluid  inside  the  control  volume  equal  the  net  mass  flow  out  of  the  control  volume 
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through  the  boundary  surface.  This  is  shown  in  equation  (2.1)  in  integral  form. 


4  [  pdV+  l  p  (V  •  h)  dS  =  0 
dt  Jn  JdQ  v  } 


(2.1) 


The  law  of  conservation  of  momentum  requires  the  time  rate  of  change  of 
momentum  due  to  unsteady  fluctuations  of  flow  properties  inside  the  control  volume 
plus  the  net  flow  of  momentum  out  of  the  control  volume  across  the  boundary  surface 
be  balanced  by  the  forces  that  act  on  the  control  volume  and  its  boundary  surface. 
This  is  shown  in  equation  (2.2)  in  integral  form. 


d_ 

dt 


J^pVdV  +  j>^pV  (V  •  n)  dS  = 


(2.2) 


Neglecting  the  body  force  terms,  the  net  force  acting  on  the  surface  of  the  control 
volume  is  given  by 

F  ~  ^  (~p  +  r)nds  (2.3) 

where  p  and  r  are  the  normal  pressure  and  shear  stress  tensor,  respectively. 

The  mass  and  momentum  integral  equations  in  conservation  form  can  be 
written  collectively  in  integral  form  as 

4fvdV+I  F  (U)  •  fids  =  0  (2.4) 

dt  Jn  Jan 

where  F  =  (Fj  — Fy)2  +  (G/“Gy)j-|-(Hj  — Hy)&.  In  differential  form  in  Cartesian 
coordinates  the  equations  can  be  written  as: 


dU  dFj  dGj  dHj  __  5Fv  dGv  dUv  ,  , 

dt  +  dx  +  dy  +  8z  dx  +  dy  +  dz  [  } 

The  dimensional  state  vector  U,  the  inviscid  flux  vectors  F/,G/,  Hj,  and  the  vis¬ 
cous  flux  vectors  Fy,Gy,  Hy,  are  expressed  in  terms  of  the  dimensional  primitive 
variables,  density  />,  pressure  p  and  x,y,z  velocities,  u ,  v,  and  w.  The  first  entry 
in  the  state  vectors  corresponds  to  the  conservation  of  mass  equation,  while  the 
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second,  third,  and  fourth  terms  correspond  to  the  x,y  and  2  momentum  equations, 
respectively. 
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2.2  Nondimensionalization 


(2.6) 

(2.7) 

(2.8) 

(2.9) 

(2.10) 

(2.11) 

(2.12) 


The  laminar  governing  equations  (2.5)  are  normalized  by  a  characteristic  length, 
Z,  freestream  density  p <*>,  and  freestream  speed,  U0 G.  Terms  in  the  equations  are 
rearranged  into  inviscid,  viscous,  and  source  terms.  The  non-dimensional  Navier- 
Stokes  equations  in  three  dimensions  can  be  written  in  Cartesian  coordinates  in 
differential  form: 


dV  &Fi  dGi  dtti_dF v  dGy  dHy 

dt  +  dx  +  dy  +  dz  dx  +  dy  +  dz  + 


(2.13) 


The  state  vector  U,  the  convective  flux  vectors  Fi,Gi,Hj,  the  viscous  flux  vectors 
Fv,  Gy,  Hy,  and  the  source  term  S  are  expressed  in  terms  of  the  non-dimensional 
primitive  variables,  namely,  density  />,  x,y,z-  velocities  and  pressure  p. 
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(2.19) 

(2.20) 


sT  =  (  0  -fl  -% 


(2.21) 

(2.22) 


The  Reynolds  number  is  defined  as 

Re=PoJJoJL 

Moo 


(2.23) 


where  is  the  freestream  density,  is  the  freestream  viscosity,  U0 0  is  the 
freestream  velocity,  and  L  is  a  characteristic  length. 


2-3  Semi-Discrete  Form  of  Governing  Equations 

An  explicit /implicit  marching  scheme  is  adopted  for  the  integration  in  time  of  equa¬ 
tions  (2.13).  The  velocity  values  are  treated  explicitly,  while  the  pressure  values 
are  treated  implicitly  in  the  momentum  equations.  The  velocity  values  are  marched 
in  time  with  a  forward  Euler  scheme  [50].  The  continuity  equation  is  formulated 
implicitly  with  the  velocity  values  considered  at  time  level  (n  +  1).  Specifically, 
the  corresponding  semi-discrete  system  is  written  as  follows,  where  the  superscripts 
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denote  the  time  levels. 
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(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 

(2.29) 

(2.30) 

(2.31) 


2.4  Pressure  Correction  Method 


Equations  (2.13)  through  (2.31)  cannot  be  solved  directly  due  to  the  implicit  treat¬ 
ment  of  the  pressure  term.  An  auxiliary  velocity  vector  U'T  =  (0,  v! ,  v9 ,  in')  is  intro¬ 
duced,  which  can  be  written  in  the  same  nondimensional  state  vector  formulation 
used  previously: 


d\J'  dFi  dGi  5Hi_5Fv  dGy  dHy 
dt  +  dx  +  dy  +  dz  dx  ^  dy  ^  dz 
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In  the  auxiliary  velocity  state  vector  equation  (2.32),  the  pressure  term  does 
not  exist,  and  U'  can  be  obtained  directly  since  it  contains  terms  only  at  the  known 
time  level  (n).  However,  the  solution  U'  does  not  satisfy  the  continuity  equation. 
Subtracting  (2.32)  and  from  the  momentum  equations  (2.13  -  2.31)  rewriting  in 
primitive  variables,  it  is  obtained: 


^+i)  -  u'  = 


At 


(2.33) 


Introducing  a  scalar  potential  <£,  such  that 


_  u’  =  -V0, 


(2.34) 


the  following  equation  for  pressure  can  be  obtained. 

P{n+1)  =  ^4>  (2.35) 

Finally,  taking  the  divergence  of  each  side  of  equation  (2.34)  and  considering  the 
continuity  equation  (2.24)  through  (2.26),  the  following  pressure  correction  Poisson 
equation  is  obtained. 

V2<t>  =  V  .  t?  (2.36) 

Using  the  <f>  values  obtained  by  equation  (2.36),  we  can  correct  the  velocity  and 
pressure  fields  using  equations  (2.34)  and  (2.35)  as  follows: 

«(n+1)  =  t?  -  V<£  (2.37) 

P{n+1)  =  ^  (2.38) 

The  present  solution  procedure  is  similar  to  the  pressure  correction  method 
in  [145].  The  overall  solution  procedure  corresponding  to  marching  by  one-time-step 
is  summarized  as  follows: 
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1.  Calculate  the  auxiliary  velocity  vector  uf  from  (2.32)  using  values. 

2.  Solve  the  pressure  correction  Poisson  equation  (2.36)  and  obtain  the  <f>  values. 

3.  Correct  the  pressure  and  velocity  at  the  n+1  time  step  using  equations  (2.37) 
and  (  2.38). 

4.  Advance  to  the  next  time  step. 
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Chapter  3 


Adaptive  Hybrid 
(Prismatic/Tetrahedral)  Grid 


The  capability  of  the  present  code  to  handle  complex  geometries  is  expanded  with 
the  use  of  a  hybrid  grid.  The  hybrid  grid  is  composed  of  semi-unstructured  prismatic 
layered  cells  near  the  surface,  and  unstructured  tetrahedra  cells  to  fill  the  remaining 
portion  of  the  domain.  Use  of  the  hybrid  grid  also  takes  advantage  of  equation 
adaptation.  The  Navier-Stokes  equations  are  solved  in  the  prismatic  region,  while 
only  the  Euler  equations  are  solved  in  the  tetrahedral  region. 


3.1  Hybrid  Mesh 

Simulation  of  flows  around  three-dimensional  bodies  is  a  major  issue  in  computa¬ 
tional  fluid  mechanics.  Geometric  and  flow-field  complexity  combine  to  make  three- 
dimensional  computations  a  pacing  item.  The  generation  of  a  body-conforming  grid 
has  proven  to  be  a  difficult  task  [11].  Success  of  structured  grid  generation  may 
be  extremely  dependent  on  geometry  and  operator  proficiency.  Block-structured 
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schemes  exist  which,  based  on  extensive  riser  input,  break  the  computational  do¬ 
main  into  a  number  of  blocks  within  which  hexahedra  are  constructed.  This  method 
requires  special  treatment  to  ensure  continuity  of  the  mesh  across  the  boundaries 
of  neighboring  blocks.  These  aspects  point  away  from  structured  grid  approaches 
toward  methods  which  require  only  suitable  definitions  of  surface  and  outer  bound¬ 
aries. 

A  radical  alternative  to  structured  meshes  is  to  use  tetrahedra.  Tetrahedral 
grids  provide  flexibility  easily  in  three-dimensional  topologies  compared  to  the  hexa- 
hedral  meshes  [76,  101,  144].  This  does  not  come  without  a  price;  unstructured  grids 
require  a  great  deal  more  memory  than  their  structured  counterparts.  They  employ 
pointers  to  provide  connectivity  information  between  cells,  faces,  edges,  and  nodes. 
Additionally,  approximately  five  to  six  times  more  tetrahedra  than  hexahedra  are 
required  to  fill  a  given  region  with  a  fixed  number  of  nodes.  This  vast  increase  in 
the  number  of  required  cells  leads  directly  to  impractical  memory  requirements  for 
three-dimensional  viscous  flow  simulations.  Furthermore,  resolution  of  the  strong 
directional  flow  gradients  encountered  in  viscous  flows  requires  very  thin  grid  ele¬ 
ments.  It  is  very  expensive  to  generate  tetrahedral  cells  with  high  aspect  ratios  to 
resolve  such  gradients. 

One  solution  to  the  dilemma  between  hexahedra  and  tetrahedra  is  to  use 
a  semi-unstructured  grid  made  of  prisms.  Prismatic  cells  are  composed  of  trian¬ 
gular  faces  in  the  lateral  (body-surface)  directions  and  quadrilateral  faces  in  the 
normal  direction  as  shown  in  Figure  3.1.  Therefore,  prisms  can  provide  the  geomet¬ 
ric  flexibility  of  unstructured  grids  as  well  as  the  orthogonality  and  high  aspect  ratio 
qualities  of  structured  grids.  Results  have  been  obtained  using  prismatic  grids  that 
reveal  their  suitability  for  resolving  viscous  flow  phenomena  [63,  66,  90,  98].  The 
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prismatic  grid  requires  a  set  of  pointers  to  define  their  base  triangular  mesh  combined 
with  a  single  index  for  each  prism  belonging  to  the  same  stack  [98].  This  reduces 
the  memory  storage  of  required  pointers  to  slightly  more  than  a  two-dimensional 
unstructured  solver. 

Areas  between  different  prismatic  layers  covering  the  surfaces  of  the  domain 
can  be  quite  irregular.  Furthermore,  relevant  flow  features  do  not  usually  exhibit  the 
strong  directionality  that  are  characteristic  of  viscous  stresses.  Tetrahedral  elements 
appear  to  be  appropriate  for  these  irregularly  shaped  regions.  Their  triangular 
faces  can  match  the  corresponding  triangular  faces  of  the  prisms.  Tetrahedral  cells 
are  created  from  these  triangular  faces  as  shown  in  Figure  3.2.  The  present  work 
employs  two  families  of  grid  elements:  prismatic  grid  cells  for  the  viscous  region 
and  tetrahedral  grid  cells  elsewhere.  In  general,  the  extent  of  the  prismatic  region 
is  specified  by  the  user  when  the  grid  is  created,  and  as  a  minimum  should  enclose 
the  entire  boundary  layer. 
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3.2  Equation  Adaptation  for  Hybrid  Mesh 


The  Navier- Stokes  equations  are  the  general  governing  equations  for  fluid  flows. 
However,  often  flows  can  be  adequately  described  with  partial  forms  of  the  Navier- 
Stokes  equations.  For  example,  inviscid  flows  can  be  adequately  computed  with  the 
Euler  equations.  The  majority  of  flows  consist  of  two  regions,  a  region  near  the  sur¬ 
face  where  the  magnitude  of  the  viscous  terms  is  appreciable,  and  an  outer  inviscid 
region  where  viscosity  can  be  neglected.  Equation  adaptation  takes  advantage  of 
this  property  and  solves  the  Navier-Stokes  equations  in  this  inner  region,  while  the 
simpler  Euler  equations  are  employed  elsewhere. 

Evaluation  of  the  viscous  terms  of  the  Navier-Stokes  equations  is  usually 
the  most  expensive  portion  of  the  computer  code.  For  incompressible  codes,  the 
solution  of  the  Poisson  equation  is  the  most  expensive  portion  of  the  code,  but  the 
resources  needed  for  the  viscous  terms  calculation  are  quite  significant.  Omitting 
the  viscous  terms  for  a  large  portion  of  the  flow  field  is  therefore  advantageous  in 
terms  of  computer  memory  storage  and  CPU  time. 

3.3  Spatial  Grid  Adaptation 

Spatial  grid  adaptation  adjusts  the  grid  by  altering  the  local  distribution  of  points. 
There  are  two  ways  of  adapting  the  grid.  The  redistribution  method  moves  the 
existing  nodes  to  cluster  the  points  in  regions  which  require  further  resolution.  The 
second  method,  embedding,  adds  additional  nodes  in  the  regions  which  require  fur¬ 
ther  resolution,  by  dividing  the  cells  in  the  region.  The  present  work  employs 
adaptive  local  embedding  of  an  initial  mesh  according  to  a  method  developed  in 

[98]. 
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3.3.1  Detection  of  Cells  to  be  Adapted 


The  success  of  spatial  grid  adaptation  is  dependent  upon  accurately  selecting  the 
regions  which  require  finer  resolution.  Regions  are  usually  selected  on  the  basis  of 
certain  flow  features.  There  are  several  flow  parameters  which  can  be  used  in  the 
flow  feature  detection  algorithm  [62,  65].  For  the  present  work,  velocity  gradients 
and  velocity  differences  will  be  the  flow  feature  parameters  used  to  determine  where 
local  embedding  should  occur.  Another  variable  in  local  adaptation  is  the  choice  of 
threshold  values  of  the  detection  parameters  for  embedding.  Regional  adaptation 
can  also  be  used  to  embed  certain  regions  of  the  mesh  based  on  x,y,z  coordinates. 

The  flow  detection  algorithm  is  edge  based.  The  variation  in  the  flow  param¬ 
eters  across  the  edges  is  monitored  during  the  adaptive  computation.  The  mean  and 
standard  deviation  of  the  velocity  gradients  and  velocity  differences  across  the  edges 
are  calculated  for  the  entire  flow  field.  If  the  value  of  the  detection  parameter,  ip, 
at  a  particular  edge  is  greater  than  the  corresponding  average  value  plus  a  fraction 
of  its  standard  deviation,  the  edge  is  flagged  for  embedding: 

V5 threshold.  =  (pmean  “h  Psd  (3.1) 

This  factor,  a,  is  called  the  threshold  parameter  and  is  determined  by  numerical 
experimentation.  The  smaller  the  threshold  parameter,  the  higher  the  percentage 
of  flagged  edges  and  more  refinement  will  occur.  If  the  threshold  parameter  is  small 
enough,  global  refinement  will  occur. 

If  two  edges  of  a  face  are  flagged  for  refinement,  the  third  edge  is  also  flagged. 
The  method  of  extended  embedding  may  also  be  used.  This  technique  introduces  ad¬ 
ditional  rows  of  embedded  faces  at  the  borders  of  the  regions  flagged  for  division.  For 
example,  if  one  edge  of  a  face  has  been  flagged  for  embedding,  extended  embedding 
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will  flag  the  remaining  two  edges  to  be  refined.  In  effect,  this  extends  the  embedded 
region  by  one  “row”  of  faces.  Extended  embedding  can  be  repeated  multiple  times 
to  flag  additional  rows  of  cells. 

The  cell- detection  procedure  is  summarized  as  follows: 

1.  Start  the  flow  simulation  on  a  relatively  coarse  grid  and  run  the  calculations 
until  the  basic  flow  features  have  developed. 

2.  Using  the  coarse  grid  and  solution  data,  loop  through  all  grid  edges  and  cal¬ 
culate  the  velocity  gradients  and  velocity  differences. 

3.  Calculate  the  mean  and  standard  deviation  of  the  distributions  of  the  velocity 
gradients  and  velocity  differences,  determine  if  the  value  of  the  velocity  gradi¬ 
ents  or  velocity  differences  of  an  edge  exceeds  the  threshold  value,  if  so  mark 
the  edge. 

4.  Apply  extended  embedding  if  desired. 

3.3.2  Local  Surface  Grid  Embedding  for  Prismatic  Mesh 

A  prismatic  mesh  is  created  from  a  triangular  mesh  distribution  on  the  surface, 
and  the  prismatic  cells  are  grown  away  from  the  surface  in  the  normal  direction 
[63,  66].  The  prismatic  grid  is  refined  by  dividing  the  triangular  surface  faces.  The 
prismatic  cells  are  then  grown  from  the  refined  surface  mesh  as  usual.  This  type 
of  adaptation  increases  the  lateral  resolution  of  the  grid,  and  does  not  destroy  the 
structure  of  the  prismatic  mesh.  A  surface  mesh  can  initially  be  relatively  coarse. 
Adaptive  embedding  subdivides  the  coarse  cells  in  the  regions  which  require  finer 
resolution.  The  subdivision  principle  is  illustrated  in  Figure  3.3.  A  triangle  can  be 
subdivided  into  two  triangles  (binary  division)  or  it  can  be  divided  into  four  triangles 
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EJ  QUAD  DIVISION 

Figure  3.3:  Local  surface  grid  embedding  for  prismatic  cells:  shaded  regions  illus¬ 
trate  binary  and  quad  triangular  face  division. 

(quad  division).  When  a  triangular  face  is  flagged  for  embedding,  it  undergoes  quad 
division.  If  the  neighbor  faces  are  not  flagged  for  refinement,  the  quad  division 
leaves  three  hanging  nodes  on  the  original  face  which  are  eliminated  by  dividing  the 
neighboring  triangular  face  into  two  triangles  (binary  division). 

The  division  of  cells  produces  new  nodes  into  the  mesh.  The  location  of 
these  new  nodes  are  interpolated  from  the  two  nodes  of  the  edge  that  is  divided.  If 
the  two  nodes  be  on  a  curved  surface,  the  interpolated  new  node  will  usually  not 
be  on  the  surface.  After  the  surface  grid  is  embedded,  the  new  nodes  are  moved  to 
a  location  which  is  consistent  with  the  curved  surface.  Interpolated  values  of  the 
velocities,  pressure,  and  pressure  correction  (cf>)  parameters  are  also  assigned  to  the 
new  nodes.  Local  embedding  can  be  repeatedly  appbed  to  obtain  multiple  levels  of 
adaptation.  This  results  in  successively  finer  locaby  embedded  grids  until  a  region 
is  adequately  resolved. 
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3.3.3  Tetrahedral  Embedding 


Tetrahedral  cells  constitute  the  portion  of  the  grid  where  inviscid  flow  features  are 
dominant.  These  features  do  not  exhibit  the  directionality  that  is  generally  prevalent 
in  viscous  stresses.  Hence,  the  tetrahedra  are  refined  isotropically,  employing  binary , 
quadtree  and  octree  divisions  as  shown  in  Figure  3.4.  As  the  feature  detector  flags 
only  the  edges  for  division,  there  will  arise  situations  where  the  flagged  edges  of 
a  tetrahedron  may  not  allow  one  of  these  three  types  of  division.  To  avoid  this 
problem  and  to  simplify  the  refinement  algorithm,  such  cases  are  constrained  by 
flagging  additional  edges  in  the  tetrahedron  so  that  the  cell  conforms  to  one  of 
the  refinement  types  shown  in  Figure  3.4.  The  tetrahedral  cells  created  by  1:8 
octree  division  have  an  aspect  ratio  that  is  comparable  to  that  of  the  parent  cell 
whereas  the  1:2  binary  and  1:4  quadtree  divisions  result  in  skewed  cells.  To  avoid 
excessive  grid  skewness,  repeated  binary  and  quadtree  divisions  of  tetrahedra  that 
originated  from  such  a  refinement  in  a  previous  pass  is  avoided.  Furthermore,  to 
avoid  sudden  changes  in  grid  size,  the  grid  refinement  algorithm  also  limits  the 
maximum  difference  in  embedding  level  between  neighboring  cells  to  two.  The 
hybrid  adaptation  is  completed  by  splicing  the  refined  prismatic  and  tetrahedral 
grids  together  at  the  interface.  The  feature  detectors  that  flag  the  edges  in  the 
prismatic  and  tetrahedral  regions  function  independently.  The  following  situations 
may  occur: 

1.  An  edge  on  the  wall  surface  may  be  flagged  for  refinement  but  its  counterpart 
in  the  tetrahedral  region  may  not  have  been  flagged  by  the  tetrahedral  feature 
detector. 
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Figure  3.4:  Refinement  strategies  for  tetrahedron:  (a)  Binary  (1:2),  (b)  Quadtree 
(1:4)  and,  (c)  Octree  (1:8)  divisions. 

2.  A  tetrahedral  edge  may  be  flagged  for  refinement,  but  its  footprint  on  the  wall 
may  not  have  been  flagged  by  the  prismatic  feature  detector. 

In  the  above  cases,  mid-edge  nodes  (hanging)  nodes  will  appear  which  would  then 
require  special  treatment  by  the  solver.  Hanging  nodes  are  avoided  by  flagging  both 
cells  of  a  prism-tetrahedron  pair  at  the  interface,  if  at  least  one  cell  is  flagged  for 


division. 


Chapter  4 


Integration  Scheme 


Spatial  integration  begins  by  conceptually  constructing  around  each  node,  N,  a  dual 
cell.  By  integrating  over  the  dual  cell,  most  of  the  integration  can  be  completed  us¬ 
ing  an  efficient  edge-based  method.  Finite  volume  discretization  of  the  momentum 
equation  is  presented  in  detail  as  is  finite  element  discretization  of  the  pressure 
Poisson  equation.  Initial  and  boundary  conditions  for  the  hybrid  grid,  artificial  dis¬ 
sipation  methodology,  time  step  calculation  and  the  solver’s  computer  requirements 
are  also  presented. 

4.1  Node-Centered  Dual  Cell 

The  present  3-D  numerical  method  uses  a  non-staggered  cell-vertex  scheme;  all 
the  primitive  variables  are  defined  at  the  cell  corners  (nodes).  The  cell  defined  by 
these  corner  nodes  is  called  the  basic  cell.  An  additional  cell  can  be  conceptually 
defined,  and  this  cell  is  called  a  node-centered  dual  cell.  The  node-centered  dual 
cell  represents  the  control  volume  over  which  the  integral  averages  of  equation  (2.4) 
are  evaluated.  The  two-dimensional  analogy  of  defining  node-centered  dual  cells 
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Figure  4.1:  2-D  analogies  for  node-centered  dual  mesh  in  the  (a)  Prismatic  region 
(b)  Tetrahedral-Prismatic  interface  and  (c)  Tetrahedral  region. 


for  different  situations  in  a  triangular-quadrilateral  hybrid  mesh  is  illustrated  in 
Figure  4.1.  Node-centered  dual  cells  are  defined  by  connecting  the  mid-points  of  the 
edges  and  centroids  of  the  triangular  and/or  quadrilateral  faces  that  share  the  node. 
Node-centered  dual  cells  for  a  three-dimensional  hybrid  grid  are  constructed  along 
similar  lines  using  the  centroids  of  faces  and  cells  that  each  node  is  associated. 

The  integral  in  equation  (2.4)  is  written  in  discrete  form  as 

=  <«> 

where  the  summation  /  is  over  all  the  discrete  faces  of  the  node-centered  dual 
mesh  that  constitute  dfijv*  The  summation  in  equation  (4.1)  can  be  alternatively 
computed  on  an  edge-wise  basis  [136]  as 

(f)K  =  ^EC =).-*.*.  <«) 

where  the  summation  e  is  over  all  the  edges  that  share  the  node  N .  The  term  S\, 
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Figure  4.2:  Node-centered  dual  cell  surfaces  attached  to  an  edge. 

represents  the  node-centered  dual  edge  area  associated  with  each  edge  and  ne  is  the 
normal  vector  of  the  node-centered  dual  edge  area  Se.  The  areas  Se  are  computed 
using  the  node-centered  dual  mesh  construction  of  Figure  4.1,  by  accumulating  the 
areas  of  each  node-centered  dual  mesh  face  that  shares  the  edge  e.  An  example  of 
the  node-centered  dual  edge  area  associated  with  a  single  edge  is  shown  in  Figure 
4.2.  The  number  of  faces  connected  to  an  edge  depends  on  the  number  of  edges  that 
shares  that  node.  Thus,  the  flux  summation  over  the  node-centered  dual  cell  surfaces 
is  equivalent  to  the  summation  over  the  edges  of  the  grid  and  the  fluxes  through 
the  node-centered  dual  mesh  faces  associated  with  each  edge.  Details  of  edge-based 
operations  are  described  in  [136].  The  finite  volume  scheme  then  proceeds  by 
computing  <5Us  at  the  nodes  by  a  global  sweep  over  the  edges  and  is  thus  transparent 
to  whether  a  node  lies  in  the  tetrahedral  region,  prismatic  region  or  at  the  interface. 


30 


4.2  Finite  Volume  Discretization  of  the  Momentum  Equa¬ 


tions 


The  Navier- Stokes  equations  of  incompressible  viscous  flow  are  discretized  using  the 
finite  volume  approach  with  a  non-staggered  grid  [61].  The  equations  are  given  in 
integral  form  for  a  bounded  three-dimensional  domain  Cl  as  follows: 


dGi 
dy  + 


dGv 

dy 


dfi 


(4.3) 


The  state,  flux  and  source  vectors  have  been  previously  defined  in  equations  (2.24) 
through  (2.31). 


4.2.1  Inviscid  Terms 

In  the  finite  volume  approach,  the  volume  integral  containing  the  spatial  derivatives 
in  equation  (4.3)  is  transformed  to  a  surface  integral  using  the  divergence  theorem. 
Considering  the  convective  flux  vector  terms,  we  have: 

hi  (I7 +  ^ + it) = hL{FlUx+GlUy  +  Hin-)dS  (44) 

Convective  fluxes  are  calculated  at  each  node,  so  the  contributions  from  each  face  of 
the  node-centered  dual  cell  is  accumulated.  The  surface  integral  on  the  right  hand 
side  of  equation  (4.4)  is  rewritten  in  discrete  form  on  an  edge- wise  basis, 

E  (FI5*  +  GI^  +  HiS2)e  (4.5) 

€ 

where  the  summation  is  over  all  the  edges,  f ln  is  the  volume  of  the  node-centered  dual 
cell,  and  Sx,  Sy,  Sz  are  the  node-centered  dual  cell  face  area  projections  associated 
with  each  edge  e  in  the  corresponding  coordinate  directions.  The  areas  Se  are 
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computed  using  the  dual  mesh  construction  of  Figure  4.1,  by  accumulating  the 
areas  of  each  node-centered  dual  cell  face  that  shares  the  edge  e.  The  flux  vectors 
are  needed  at  the  center  of  the  edge  e  and  their  values  are  obtained  by  averaging 
the  values  at  the  two  nodes  of  each  edge. 


4.2.2  Viscous  Terms 


The  volume  integral  containing  the  spatial  derivatives  in  equation  (4.3)  is  trans¬ 
formed  to  a  surface  integral  using  the  divergence  theorem.  Considering  the  viscous 


flux  vector  terms,  we  have: 


dx  dy 


=  7 r  /  (Fynj;  +  Gv%  +  (4-6) 


The  flux  vectors  Fy,Gv,  and  Hv  contain  the  gradients  of  the  velocities  that  need 
to  be  calculated  first  at  the  edge  centers.  For  this  edge  center  calculation,  another 
conceptual  dual  cell  is  used  [18].  The  two-dimensional  edge-centered  dual  cell  is 
shown  in  Figure  4.3.  The  edge-centered  dual  cell  is  composed  of  all  the  cells  which 
share  the  edge.  For  an  unstructured  grid,  the  number  of  cells,  m,  which  share  an 
edge  will  vary. 

Equation  (4.7)  shows  the  calculation  for  the  term. 


^  =  2-  [  unxdS 

fie  JU  OX  fie  JdQ 


i  7TL 

cT 


This  summation  is  over  the  faces  of  the  edge-centered  dual  cell,  which  correspond  to 
selected  faces  of  the  basic  cells.  fie  is  the  volume  associated  with  the  edge-centered 
dual  cell,  and  SxySy,  and  Sz  are  the  area  projections  of  the  edge-centered  dual 
cell  face  /  in  the  corresponding  coordinate  directions.  Velocities  are  required  at  the 
center  of  the  edge-centered  dual  face  /  and  their  values  are  obtained  by  averaging  the 
values  at  the  four  nodes  of  the  quadrilateral  face  or  three  nodes  of  the  triangular  face. 
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Figure  4.3:  2-D  analogies  for  edge-centered  dual  mesh  in  the  (a)  Prismatic  region 
(b)  Tetrahedral-Prismatic  interface  and  (c)  Tetrahedral  region. 

For  the  hybrid  grid,  when  inviscid  tetrahedra  are  used,  the  summation  is  over  the 
edge-centered  dual  cell  faces  in  the  prismatic  region.  When  viscous  tetrahedra  are 
used,  the  summation  is  over  the  edge-centered  dual  cell  faces  in  the  entire  domain. 

The  second  derivatives  (divergence  of  the  velocity  gradients)  are  calculated  at 
the  nodes.  The  right  hand  side  of  equation  (4.6)  is  discretized  over  the  node-centered 
dual  cells  and  equation  (4.8)  is  obtained.  The  integration  is  completed  using  the 
node-centered  dual  mesh  and  edge-based  methodology  described  in  the  previous 
section  on  the  inviscid  terms.  For  the  hybrid  grid,  when  inviscid  tetrahedra  are 
used,  the  summation  is  over  the  prismatic  edges.  When  viscous  tetrahedra  are  used, 
the  summation  is  over  all  the  edges. 

E  (FvS*  +  GvSy  +  Hv5,)e  (4.8) 

^71  P 
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4.3  Finite  Element  Discretization  of  Poisson  Equation 


Finite  element  discretization  of  the  Poisson  equation  is  used  because  the  resulting 
coefficient  matrix  is  reasonably  diagonally  dominant.  Tetrahedral  finite  elements  are 
chosen  since  these  elements  have  attractive  properties  and  the  coefficient  matrix  is 
simpler  to  construct  than  any  other  element.  The  hybrid  mesh  is  already  composed 
of  a  majority  of  tetrahedral  cells,  and  only  requires  conceptual  division  of  the  prism 
cells  into  tetrahedra  cells  for  the  formation  of  the  Poisson  equation  matrix.  This 
prismatic  division  procedure  is  described  in  Appendix  A. 

4.3.1  Isoparametric  Transformation 

The  pressure  correction  equation  (2.36)  is  discretized  using  the  Galerkin  finite- 
element  approach  [12].  The  scheme  is  compact  with  all  operations  being  restricted 
to  within  each  grid  cell.  Bilinear  isoparametric  tetrahedral  elements  are  employed 
[39,  148].  For  tetrahedral  cells,  the  values  of  </>  and  k  are  defined  in  each  element 
using  the  following  finite-element  formulation  : 

4  4 

cf>  =  J2^i,  k  =  '£N1k1  (4.9) 

4=1  4  =  1 

where  (f>{  are  nodal  values  of  </>,  K{  are  arbitrary  constants  of  the  test  function  k,  and 
N{  is  the  interpolating  shape  function  associated  with  the  z-th  node  in  x,y,z  space. 
The  x,  y,  and  z  coordinates  can  also  be  defined  in  terms  of  the  shape  functions: 

4  4  4 

x  =  NiXi,  y  =  ^2  Nm,  z  =  J2  Niz*  (4.10) 

4=1  4  =  1  4=1 

The  master  tetrahedral  element  is  shown  in  Figure  4.4  and  the  bilinear  shape 
functions  for  the  master  tetrahedral  element  in  sPace  are  defined  as 
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Figure  4.4:  Master  tetrahedral  element  in  £,?7,C  space. 


Nl  =  (4.11) 

N2  =  t 

iV3  =  7? 

n4  =  c 


Relations  are  needed  to  relate  the  master  shape  functions  in  £,//,£  space  to  the 
corresponding  shape  functions  in  x,  y,  z  space.  The  transformation  needed  to  supply 
these  relations  [148]  is  shown  in  equation  (4.12). 


dN 
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dz 

(4.12) 


The  matrix  in  equation  (4.12)  is  the  Jacobian  of  the  transformation.  The 
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derivatives  of  the  shape  functions  in  x,y,  and  z  space  can  be  found  by  inverting 
equation  (4.12). 

All  the  terms  of  the  Jacobian  matrix,  can  be  calculated  directly  since  the 
relationship  between  physical  coordinates  and  shape  functions  are  known  from  equa¬ 
tions  (4.10)  through  (4.11).  Shape  function  derivatives  can  be  found  easily  since  the 
shape  functions  are  defined  in  £  ,7/ ;  C  space.  The  Jacobian  matrix  is  a  3x3  matrix 
and  can  be  easily  inverted  in  the  numerical  code. 


4.3.2  Galerkin  Method 

The  Galerkin  method  is  applied  to  the  left  hand  side  of  equation  (2.36)  to  obtain 
the  matrix  coefficients  of  the  Poisson  equation.  Rewriting  the  equation,  we  have: 

-A  <f>  =  f(x,y,z)  (4.13) 


where  f(x,y,z )  is  the  divergence  of  the  velocity.  Multiplying  both  sides  by  the  test 
function  k,  and  integrating  over  the  domain  Cl  we  obtain: 

/  (-kA</>  —  nf)  dCl  =  0  (4.14) 

Jn 

Applying  the  following  vector  identity  to  equation  (4.14) 

-  kA <f>  =  -V  •  (kV<£)  +  V<£  •  Vk  (4.15) 


we  obtain 


[  [V<£  •  Vk  -  V  •  (kV</>)  -  nf]  dCl  =  0 

Jn 


(4.16) 


Using  the  divergence  theorem  on  the  second  term  on  the  left  hand  side  of  equation 
(4.16),  we  obtain  equation  (4.17). 


/  _v  •  (KV<t>)d,Cl  =  -  [  ^ KdS 
Jn  Jan  dn 


(4.17) 


36 


This  term  will  disappear  since  the  integral  of  over  the  surfaces  of  consistently 
defined  tetrahedral  cells  is  zero.  Therefore  the  pressure  correction  equation  becomes 

/  f  KfdO,  (4.18) 

Ju  Jn 

The  coefficients  of  <f>,  determined  from  the  left  hand  side  of  the  equation,  become 
the  entries  in  the  Poisson  equation  matrix. 

The  element  matrix  system  therefore  becomes 


Dij<t>j  —  fj 


(4.19) 


where 


A;  = 


dNj  dNj 
dx  dx 


+ 


dNj  dNj 
dy  dy  + 


dNj  dNj 

dz  dz 


^  dx  dy  dz 


(4.20) 


and  the  derivatives  in  equation  (4.20)  are  found  after  inverting  the  Jacobian  in 
equation  (4.12).  The  integrand  is  constant  over  each  element  and  the  integration 
can  be  completed  by  multiplying  the  factor  in  the  parentheses  with  the  volume  of 
the  tetrahedral  cell. 


Instead  of  the  finite  element  technique,  the  right  hand  side  vector,  fj,  is 
calculated  using  the  finite  volume  technique,  similar  to  the  edge-based  method  shown 
earlier  for  the  inviscid  momentum  equation  terms. 


fj  =  '%2(uSx  +  vSy  +  wSz)e  (4.21) 

e 


4.3.3  Element  and  Global  Matrix  Construction 


We  can  construct  the  following  global  matrix  system  by  assembling  the  element 
matrix  obtained  from  equation  (4.19) 


D$  =  f 


(4.22) 
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where 


D  =  ]T  Dij  (4.23) 

e 

f  =  E/. 

e 

$T  =  [^1,^2,  ^3,‘"] 

Values  for  <f>  on  the  left  hand  side  of  equation  (4.22)  are  treated  implicitly,  which 
requires  the  solution  of  a  system  of  equations.  The  matrix  that  is  formed  after  these 
steps  is  a  symmetric  linear  matrix.  This  matrix  is  slightly  modified  to  handle  the 
implicit  boundary  conditions,  and  in  general  is  no  longer  symmetric. 

Solution  of  the  pressure  Poisson  equation  is  the  most  computationally  ex¬ 
pensive  portion  of  the  incompressible  flow  calculation  [81].  The  fraction  of  com¬ 
putational  time  consumed  for  solving  this  equation  can  be  as  high  as  80  %  of  the 
total  CPU  cost  [134].  It  is  very  important  that  the  resulting  system  of  equations 
be  solved  in  an  efficient  manner.  The  use  of  direct  solvers  is  unattractive  because 
of  the  large  storage  requirements  and  computer  effort.  The  matrix  created  by  the 
Poisson  equation  is  large  and  sparse.  The  full  matrix  is  an  n  by  n  matrix,  where  n 
is  the  total  number  of  nodes.  The  system  is  solved  by  the  Gauss-Seidel  method. 

Use  of  iterative  solvers  creates  the  question  of  when  the  iteration  procedure 
should  be  terminated.  Often  the  convergence  criterion  is  chosen  as  an  arbitrarily 
small  number.  Performance  of  the  entire  algorithm  depends  on  the  convergence 
criteria  for  the  Poisson  equation.  If  the  Poisson  iteration  is  terminated  before  suffi¬ 
cient  convergence  is  achieved,  errors  are  propagated  throughout  the  solution,  possi¬ 
bly  resulting  in  divergence  or  slow  convergence.  However,  it  is  a  waste  of  computer 
resources  to  choose  convergence  criteria  which  are  too  small. 

Different  values  of  the  Poisson  equation  convergence  criteria  were  used  to 
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determine  its  effect  on  the  solution.  Increasing  the  Poisson  convergence  criteria 
significantly  increased  the  number  of  iterations  needed  to  reduce  the  divergence  to 
an  acceptable  level.  In  the  present  work,  when  the  absolute  value  of  the  difference 
in  the  value  of  <f>  between  iterations  is  less  than  or  equal  to  5.0e-5  for  any  node  i, 
convergence  is  reached.  This  convergence  criteria  is 

<  5.0e-  5  (4.24) 

where  k  denotes  the  iteration  level. 

In  the  solution  of  the  Poisson  equation,  it  has  been  found  the  lagging  of 
boundary  values  is  detrimental  [103].  Boundary  values  must  be  updated  in  a 
manner  similar  to  the  interior  points.  Therefore,  the  present  work  uses  the  iterative 
solver  to  solve  for  all  the  grid  points  simultaneously. 

4.4  Boundary  Conditions  for  Unstructured  Grids 

Treatment  of  boundary  conditions  is  crucial  to  any  numerical  simulation.  With 
structured  meshes,  the  boundary  conditions  can  be  treated  easily  using  the  i,j,k 
indices.  In  unstructured  meshes,  treatment  of  the  boundary  conditions  are  not  as 
straightforward.  Treatment  of  surface,  far  field,  symmetry  planes  and  side  walls 
developed  in  the  present  work  is  presented. 

4.4.1  Wall 

Surface  wall  boundary  conditions  are  applied  to  the  three  velocity  components,  as 
well  as  to  the  pressure  corrections,  (f>.  At  the  wall,  the  u ,  v  and  w  components  of 
velocity  are  set  to  zero.  The  value  of  d(f)/dn  is  also  set  to  zero. 
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Figure  4.5:  Obtaining  the  number  of  the  node  adjacent  to  the  wall  node 

Since  the  wall  boundary  nodes  are  all  located  on  triangular  faces  of  prismatic 
cells,  the  boundary  condition  d(f>/dn  can  be  treated  very  simply.  Prismatic  layers 
are  grown  away  from  the  wall  surface  in  the  normal  direction  so  the  neighbor  node 
which  is  normal  to  a  specific  wall  node  is  known  immediately.  For  example,  suppose 
the  total  number  of  nodes  on  the  wall  boundary  is  ibnodtot.  The  node  immediately 
off  the  wall  node,  j,  in  the  normal  direction  is  obtained  by  adding  ibnodtot  to  the 
wall  node  number,  as  shown  in  Figure  4.5.  The  wall  boundary  is  the  shaded  region. 
In  the  process  of  solving  the  Poisson  equation  and  satisfying  dcj)/dn  =  0  at  node  j, 
we  prescribe,  ^(j)  =  </>(j+ibnodtot).  This  is  completed  for  all  the  nodes  on  the  wall 
boundary. 

4.4.2  Far  Field 

The  far  field  boundary  conditions  are  applied  to  the  three  velocity  components,  as 
well  as  to  the  pressure  corrections,  At  the  far  field ,  the  velocity  components  are 
set  equal  to  the  freestream  values,  while  the  value  of  d<f>/dn  is  set  to  zero.  The 
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treatment  of  the  far  field  boundary  conditions  of  the  prismatic  mesh  are  different 
than  the  treatment  of  boundary  conditions  of  the  hybrid  mesh. 

Prismatic  Grid 

For  the  all-prismatic  grid,  the  far  field  boundary  nodes  are  all  located  on  triangu¬ 
lar  faces  of  the  prismatic  cells,  and  the  boundary  condition  d(j)/dn  can  be  treated 
similarly  to  the  treatment  of  the  surface  wall  boundary  nodes.  In  the  process  of 
solving  the  Poisson  equation  and  satisfying  d(f>/dn  =  0  at  far  field  node  j,  <^(j)  = 
</>(j-ibnodtot)  is  prescribed.  This  is  completed  for  all  the  nodes  on  the  far  field 
boundary. 

Hybrid  Grid 

In  the  hybrid  grid,  the  tetrahedra  cells  are  grown  away  from  the  outer  prismatic 
region  to  fill  a  hexahedral  domain.  The  far  field  boundary  lies  along  the  x  =  ± 
xmax,  y  =  ymax,  and  z  =  ±  zmax  planes.  The  far  field  boundary  nodes  in  a  hybrid 
grid  are  triangular  faces  of  a  tetrahedra.  When  a  tetrahedral  cell  is  on  the  far  field 
boundary,  the  cell  center  value  of  </>,  is  extrapolated  to  the  center  of  the  triangular 
face  on  the  boundary.  This  face  center  value  is  distributed  to  the  face  nodes.  The 
nodal  distribution  uses  weighting  factors  which  are  inversely  proportional  to  the 
distance  the  face  node  is  away  from  the  face  center.  Therefore,  the  closer  the  face 
node  is  to  the  face  center,  the  greater  the  distribution  the  face  node  will  receive. 
This  procedure  is  illustrated  in  Figure  4.6,  where  the  far  field  boundary  surface  is 
the  shaded  region.  In  the  Poisson  equation  calculation,  this  methodology  is  treated 
implicitly. 

Care  must  be  taken  when  defining  the  faces  on  the  far  field  boundary  of  the 
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Figure  4.6:  Illustration  of  far  field  boundary  conditions  for  a  hybrid  mesh:  assigning 
the  cell  center  value  to  the  tetrahedral  face  center,  and  distributing  the  face  center 
value  to  the  face  nodes. 

hybrid  grid.  A  face  is  defined  as  being  on  the  far  field  boundary  if  its  three  nodes  all 
lie  on  the  same  far  field  boundary  plane .  A  face  can  have  all  three  nodes  identified 
as  a  far  field  boundary  node,  yet  the  face  defined  by  the  three  nodes  may  not  lie  on 
the  boundary.  This  is  easily  seen  for  a  tetrahedral  cell  which  lies  at  the  juncture 
of  two  far  field  planes  of  the  computational  domain,  as  shown  in  Figure  4.7.  The 
shaded  planes  lie  on  the  far  field  boundary. 

4.4.3  Symmetry  Planes 
Prismatic  Grid 

Boundary  nodes  on  a  symmetry  plane  in  the  prismatic  region  are  nodes  on  a  quadri¬ 
lateral  face.  Since  the  cells  are  prisms,  there  is  no  direct  one-to-one  correspondence 
of  nodes  normal  to  the  nodes  on  the  symmetry  plane.  The  boundary  conditions  on 
a  symmetry  plane  are  treated  in  a  methodology  similar  to  the  hybrid  far  field  nodes. 
Symmetry  boundary  calculations  are  performed  in  the  following  manner. 
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Figure  4.7:  Illustration  of  a  face  defined  by  three  far  field  nodes  which  does  not  lie 
on  the  far  field  boundary  (shaded  region). 

The  cell  center  value  of  the  x,  y,  z- velocity  or  (f>  is  determined  first  by  averaging  the 
values  from  the  six  cell  nodes.  The  cell  center  value  is  extrapolated  to  the  center 
of  the  face  on  the  symmetry  plane.  The  face  center  value  is  then  distributed  to 
the  face  nodes.  The  nodal  distribution  uses  weighting  factors  which  are  inversely 
proportional  to  the  distance  the  face  node  is  away  from  the  face  center.  Therefore, 
the  closer  the  face  node  is  to  the  face  center,  the  greater  the  distribution  the  face 
node  will  receive.  This  procedure  is  illustrated  in  Figure  4.8.  In  the  figure,  the 
symmetry  boundary  surface  is  the  shaded  region. 

In  the  case  of  the  velocities  at  a  symmetry  plane,  the  velocity  normal  to  the 
plane  is  set  to  zero.  The  values  of  the  other  velocity  components  are  treated  as 
described  above.  In  the  case  of  the  pressure  correction  on  a  symmetry  plane,  the 
values  of  (f)  at  the  nodes  are  handled  implicitly  and  the  subsequent  equation  becomes 
part  of  the  Poisson  matrix. 

Boundary  conditions  for  the  symmetry  planes  for  the  cylinder  case  were 
slightly  different  since  the  cylinder  has  a  strictly  two-dimensional  flow.  Instead  of 


Figure  4.8:  Illustration  of  symmetry  boundary  conditions:  assigning  the  cell  center 
value  to  the  quadrilateral  face  center,  and  distributing  the  face  center  value  to  the 
face  nodes. 

extrapolating  the  <j>  cell  center  value  to  the  face  center,  the  (f>  value  of  the  (third)  node 
on  the  triangular  face,  not  on  the  symmetry  edge,  is  directly  distributed  to  the  two 
nodes  on  the  edge,  using  a  weighting  factor  inversely  proportional  to  the  distance  the 
edge  node  is  away  from  the  projection  of  the  third  node  on  the  edge.  This  reduces 
the  transportation  of  <f>  values  in  the  radial  and  circumferential  direction.  For  the 
cylinder  case,  0  values  should  only  be  transferred  in  the  axial  (normal)  direction  for 
a  true  symmetry  boundary  condition. 

Hybrid  Grid 

Symmetry  boundary  conditions  for  the  tetrahedral  portion  of  the  hybrid  grid  are 
handled  identically  to  the  far  field  boundary  conditions  of  the  hybrid  grid.  The 
symmetry  plane  for  the  hybrid  grid  lies  at  y=0.0.  When  a  tetrahedral  cell  is  on  the 
symmetry  boundary,  the  cell  center  value  of  the  velocity  or  <j>,  is  extrapolated  to  the 
center  of  the  triangular  face  on  the  boundary.  This  face  center  value  is  distributed 
to  the  face  nodes.  The  nodal  distribution  uses  weighting  factors  which  are  inversely 
proportional  to  the  distance  the  face  node  is  away  from  the  face  center.  Therefore, 
the  closer  the  face  node  is  to  the  face  center,  the  greater  the  distribution  the  face 


node  will  receive.  In  the  Poisson  equation  calculation,  this  methodology  is  treated 
implicitly. 

4.5  Initial  Conditions 

Initial  conditions  for  the  flow  around  a  cylinder  and  sphere  were  zero  velocity  at  the 
surface  and  a  linear  distribution  of  velocity  to  the  outer  edge  of  the  prisms  where  the 
velocity  was  set  equal  to  the  freestream  value.  For  the  hybrid  mesh,  the  velocity  in 
the  tetrahedral  region  is  initially  the  freestream  value.  The  pressure  (oc  pressure  co¬ 
efficient)  was  assumed  to  be  a  Bernoulli  equation  type  pressure  distribution  initially, 
p  =  \u2.  Initial  conditions  for  <fi  were  pSt ,  where  6t  was  an  initial  representative 
value  of  0.0001. 

Initial  conditions  for  the  flow  over  a  flat  plate  were  zero  pressure  and  freestream 
velocity  everywhere  except  for  no  slip  conditions  at  the  wall  and  the  specified  ve¬ 
locity  at  the  inlet.  Initial  conditions  for  the  flow  inside  a  driven  cubic  cavity  were 
stagnant  flow  conditions  except  for  the  flow  at  the  moving  wall.  <f>  was  assumed  to 
be  zero  initially. 


4.6  Artificial  Dissipation 

Central  space  differencing  schemes  are  susceptible  to  oscillatory  modes  in  the  veloc¬ 
ity  field  of  high  Reynolds  number  flows.  Furthermore,  odd-even  decoupling  of  the 
solution  may  appear  in  the  pressure  field  for  this  non-staggered  type  of  mesh  that 
is  employed.  To  stabilize  the  calculations  for  high  Reynolds  number  flows,  artifi¬ 
cial  dissipation  is  often  used  explicitly  or  implicitly  [69].  For  low  Reynolds  number 
flows,  artificial  dissipation  may  not  be  required  since  the  viscous  terms  are  relatively 
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large  and  can  eliminate  the  error  terms  which  have  been  introduced.  In  the  present 
work,  a  fourth  order  smoothing  term  is  added  explicitly  to  the  momentum  equations 
in  order  to  suppress  odd-even  decoupling  of  the  solution  [57). 

The  smoothing  operator  is  cast  in  a  form  suitable  for  adaptive  unstructured 
grids.  It  is  formulated  in  such  a  manner  as  to  simulate  the  implicit  dissipation 
terms  of  the  upwinding  schemes  without  increasing  the  computational  cost  of  the 
algorithm  [104,  136).  Edge-based  operations  are  used  for  the  calculation  of  the 
artificial  dissipation  term.  Consider  an  edge  formed  by  the  nodes  L  and.  R. 

The  equation  for  the  second  order  contributions  at  node  L  is  shown  in  equa¬ 
tion: 

(*UL)2  =  £(Ufi-UL)  (4.25) 

e=l 

where  e  denotes  the  edges  sharing  node  L.  The  fourth  order  smoothing  contribution 
is  computed  in  a  similar  fashion.  Instead  of  the  first  difference  of  state  vectors  as 
used  in  equation  (4.25),  a  difference  of  the  accumulated  first  difference  over  the 
edges  sharing  a  node  is  used  for  the  fourth  order  smoothing  contribution.  The 
fourth  order  difference  is  scaled  by  the  time  step  at  node  L,  At,  and  node  volume, 
$lL.  This  scaling  factor  was  determined  by  numerical  experimentation  and  was  the 
most  effective  of  the  scaling  factors  tried. 

A  /  n 

(«UL)4  =  -a4—  £  [(«SU*)2  -  («JL)J  (4.26) 

llL  e=l 

The  coefficient  <74  is  an  empirical  parameter  that  controls  the  amount  of  fourth  order 
smoothing. 

Large  values  of  <t4w,  ct4v,  or  <74^,  stabilize  the  solution  but  destroy  the  accu¬ 
racy.  Therefore,  special  care  is  required  for  choosing  the  value  of  the  smoothing  coef¬ 
ficients.  Numerical  experiments  have  been  carried  out  to  determine  optimum  values 
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for  the  smoothing  coefficients.  The  determined  values  are  such  that  the  solution 
accuracy  is  not  affected,  while  the  odd-even  modes  are  suppressed  [61].  Additional 
dissipation  in  the  pressure  Poisson  equation  was  not  used.  Numerical  experimenta¬ 
tion  showed  the  artificial  dissipation  in  the  momentum  equation  induced  adequate 
smoothing  in  the  pressure  field. 


4.7  Time  Step  Calculation 


The  stability  limitation  for  the  model  1-D  convection  equation  ut  +  cuy  =  0  is 
^  <  1  (CFL  limitation),  while  the  corresponding  stability  restriction  for  the  1-D 
model  diffusion  equation  ut  =  vuyy  is 

In  the  present  central  space  and  forward  time  differenced  scheme,  a  combi¬ 
nation  of  the  two  limitations  is  employed.  Specifically, 


A*  = 


u 


V0h 


|  U{  Six  |  4"  l^i  Siy  |  \W{Siz  |  + 


2'volj 


(4.27) 


where  vol{  is  the  dual  volume  corresponding  to  node  i,  and  W{  are  the  velocities 

at  node  i,  Six,  S{y,  and  S{z  are  the  projected  areas  of  the  dual  faces  for  the  dual 
cell  at  node  i,  and  Re  is  the  Reynolds  number.  Lastly,  to  is  the  CFL  (Courant, 
Friedrichs,  and  Lewy)  factor  [7],  and  equal  to  a  number  around  0.1.  The  small  value 
for  v  reflects  the  property  of  projection  methods  requiring  a  restrictive  time  step 


[13]- 


4.8  Solver  Computer  Requirements 

The  incompressible  hybrid  solver  is  described  in  terms  of  computing  time  and  mem¬ 
ory  requirements.  The  CPU  time  consumed  on  a  IBM  390  workstation  is  0.0002 
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seconds/node/time  iteration.  Table  4.1  illustrates  the  percentage  of  CPU  time  which 
is  consumed  by  the  various  parts  of  the  solver.  These  percentages  were  obtained 
from  averaging  the  CPU  times  consumed  by  the  first  100  iterations  of  a  simulation 
of  flow  over  a  sphere  using  a  hybrid  mesh.  The  pressure  correction  procedure  is 
the  most  expensive  in  terms  of  CPU  time,  because  it  contains  the  solution  of  the 
Poisson  equation  matrix  using  the  Gauss- Seidel  method.  The  first  ten  time  steps 
of  the  simulation  requires  many  iterations  of  the  Gauss- Seidel  method,  however, 
most  of  the  simulation  use  one  Gauss-Seidel  iteration  per  time  step.  The  percent¬ 
age  of  time  for  the  pressure  correction  method  would  decrease  if  data  from  more 
iterations  were  averaged.  The  percentage  of  CPU  time  which  is  consumed  by  the 
pressure  correction  part  of  the  solver  is  further  examined  in  Table  4.2.  The  viscous 
subroutine  is  the  next  most  expensive  subroutine.  For  this  hybrid  mesh,  equation 
adaptation  was  used.  There  are  19740  nodes  in  the  prismatic  region  which  are  in¬ 
volved  in  the  Navier-Stokes  calculation.  There  are  another  10418  nodes  solely  in  the 
tetrahedral  region  which  are  involved  in  the  Euler  calculation  and  not  involved  in 
the  viscous  calculation.  If  the  Navier-Stokes  equations  were  solved  over  the  entire 
domain,  the  amount  of  time  consumed  by  the  viscous  calculations  would  increase 
by  30  %,  thereby  increasing  its  percentage  of  CPU  time  consumed.  It  is  clear  from 
examination  of  CPU  time,  there  is  a  definite  advantage  in  utilizing  equation  adap¬ 
tation. 

The  solver  is  written  in  a  form  which  permits  vectorization  on  vector  parallel 
machines.  Most  of  the  computationally  significant  subroutines  can  be  vectorized 
with  the  exception  of  the  Gauss-Seidel  subroutine,  resulting  in  significant  reduction 
in  CPU  time.  It  is  possible  to  force  the  vectorization  of  the  Gauss-Seidel  subrou¬ 
tine,  but  the  resulting  matrix  solver  would  probably  exhibit  convergence  properties 
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Table  4.1:  Proportions  of  computing  time  consumed  by  parts  of  the  solver. 


Subroutine 

CPU  Time  Proportion  (%) 

Convective  Terms 

7.0 

Viscous  Terms 

33.6 

Smoothing 

2.5 

Boundary  Conditions 

0.5 

Time  Step  Calculation 

0.1 

Pressure  Correction 

56.3 

between  Jacobi  and  Gauss-Seidel  techniques.  The  savings  in  CPU  time  through 
vectorizing  may  offset  the  additional  iterations  needed  from  the  resulting  mixed 
Jacobi-Gauss-Seidel  solver. 

Memory  requirements  for  the  solver  are  shown  in  Table  4.3.  The  variables 
requiring  the  most  memory  are  those  associated  with  the  Poisson  matrix  and  vectors. 
The  bytes  shown  for  the  Poisson  matrix  are  representative  numbers  for  a  prismatic 
mesh  where  nodes  have  a  maximum  of  20  neighboring  nodes.  This  is  a  generous 
number  for  a  prismatic  mesh.  However,  a  hybrid  mesh  may  result  in  four  times 
as  many  neighbors,  thus  substantially  increasing  the  amount  of  storage  needed  for 
the  Poisson  equation.  The  majority  of  Poisson  equation  variables  are  also  double 


49 


Table  4.2:  Proportions  of  computing  time  consumed  by  pressure  correction  parts  of 
the  solver. 


Subroutine 

CPU  Time  Proportion  (%) 

Divergence  of  Velocity 

4.7 

Gauss-Seidel 

88.3 

Gradient  of  Phi 

5.7 

Boundary  Conditions 

0.8 

precision,  requiring  twice  as  many  bytes  compared  to  the  other  storage  arrays.  A 
smaller,  but  still  significant,  memory  storage  “sink”  are  the  pointers  and  metrics 
needed  for  the  tetrahedra  cells  if  the  Navier-Stokes  equations  are  solved  in  this 
region.  Again,  the  advantage  to  using  equation  adaptation  is  clear,  even  when 
based  on  memory  requirements  alone. 
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Table  4.3:  Memory  requirements  of  the  solver. 


Item 

Bytes 

State  Vectors:  U  and  <5U,  per  node 

40 

Pressure  and  Time  Step,  per  node 

8 

Poisson  Matrix  and  Vectors,  per  node 

272 

Node  Metrics,  per  node 

28 

Edge  Metrics,  per  edge 

24 

Boundary  Face  Metrics,  per  boundary  face 

112 

Boundary  Edge  Metrics,  per  boundary  edge 

32 

Viscous  Tetrahedra,  per  tetrahedral  face 

52 
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Chapter  5 


Numerical  Issues 


In  the  course  of  the  development  of  the  present  numerical  method,  several  numerical 
issues  arose  which  had  significant  impact  on  the  direction  of  the  research.  These 
issues  include  the  pressure  correction  formulation,  diagonal  dominance  of  the  result¬ 
ing  pressure  Poisson  equation,  boundary  conditions,  non-staggered  versus  staggered 
grids,  discrete  divergence,  and  testing  subroutines  with  analytic  field  functions. 

5.1  Pressure  Correction  Formulations 

In  the  evolution  of  the  present  method’s  development  I,  as  well  as  others  [109], 
have  found  that  the  particular  formulation  of  the  Poisson  pressure  equation  has  a 
tremendous  impact  on  the  quality  of  the  solution. 

The  pressure  correction  methods  tested  in  the  present  work  are  variants  of 
Chorin’s  original  projection  method  which  was  based  on  the  Hodge  projection  [30]. 
Projection  methods  are  a  type  of  fractional  step  method  [14,  30,  145].  These 
methods  obtain  an  auxiliary  velocity  vector  field  by  disregarding  the  solenoidal 
nature  of  the  velocity,  and  then  project  this  auxiliary  vector  field  onto  the  space  of 
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divergence  free  fields  to  obtain  the  divergence  free  velocity  [13,  30].  The  numerical 
form  of  the  projection  method  computes  an  auxiliary  vector  field  and  applies  some 
type  of  discrete  projection.  The  projection  recovers  the  divergence  free  velocity  by 
splitting  the  velocity  field  into  two  parts  that  are  divergence  free  and  curl  free  [108]. 
The  curl  free  part  is  written  as  the  gradient  of  a  potential.  The  auxiliary  vector 
field  is  a  linear  combination  of  the  velocity  at  the  new  time  level  and  the  gradient 
of  the  potential. 

u'  =  +  V0 


The  first  method  tried  was  the  Simplified  Marker  and  Cell  (SMAC)  [4,  5,  71] 
method,  a  modification  of  the  Marker  and  Cell  (MAC)  [50]  method.  The  SMAC 
solution  procedure  follows  the  explicit/implicit  marching  scheme  in  [71,  91].  The 
overall  solution  procedure  is  very  similar  to  the  method  chosen  in  the  present  work 
except  in  the  choice  of  the  auxiliary  velocity  vector.  The  current  auxiliary  velocity 
vector  was  defined  in  equation  (2.32)  and  is  repeated  here  for  clarity. 

^  +  ^  +  +  «Hi  =  MV  dGv  3Hv  (  * 

dt  dx  dy  dz  dx  dy  dz 


In  the  SMAC  method,  the  auxiliary  velocity  vector  equation  contains  the  pressure 
at  the  n  time  level,  equation  (5.2). 


dW  dFi  dGi  dHi  dlW  dGv  dHv 
dt  +  dx  dy  +  dz  dx  ^  dy  dz  + 


(5.2) 


Similarly,  the  auxiliary  velocity  state  vector  can  be  obtained  directly  since 
it  contains  terms  only  at  the  known  n  time  level,  and  it  does  not  satisfy  the  con¬ 
tinuity  equation.  Subtracting  the  auxiliary  velocity  vector  equation  (5.2)  from  the 
momentum  equation  (2.13)  through  (2.31),  and  rewriting  in  primitive  variables,  it 
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is  obtained  : 


tS<n+1>  -  u'  =  -V  Un+1>  -  pH  At  (5.3) 

Introducing  a  similar  scalar  potential  0,  such  that 

U<n+1)  _  u'  =  -V<f>,  (5.4) 

the  following  equation  for  pressure  can  be  obtained: 

p(n+D_p(n)  =  ^  (5.5) 

Finally,  taking  the  divergence  of  each  side  of  equation  (5.4)  and  considering  the 
continuity  equation  (2.24)  through  (2.26),  the  following  pressure  correction  Poisson 
equation  is  obtained  for  the  SMAC  formulation. 

A<j>  =  V-t?  (5.6) 

Using  the  (j)  values  obtained  by  equation  (5.6),  we  can  correct  the  velocity  and 
pressure  fields  using  equations  (5.4)  and  (5.5),  derived  from  the  SMAC  method,  as 
follows: 

rfn+l)  _  _  V(£  (5<7) 

p(«+D=J,(n  )  +  _L^  (5.8) 

The  major  difference  between  the  present  method  and  the  SMAC  method 
is  the  pressure  correction  equation,  equations  (2.38)  and  (5.8).  Equation  (2.38), 
repeated  here  for  convenience, 

P(n+1)  =  ^  (5.9) 

obtains  the  corrected  pressure  field  from  0,  while  equation  (5.8)  yields  the  cor¬ 
rected  pressure  field  from  the  previous  pressure  field  and  the  correction.  The  SMAC 
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Figure  5.1:  Wall  pressure  coefficient  on  cylinder  for  Re=40. 
o  Experiment  by  Grove  [46] 

- SMAC  method 

- Present  work 

method  accumulates  the  pressure  changes.  As  the  iterative  accumulation  took  place 
using  the  SMAC  method,  an  oscillatory  pressure  field  developed.  The  wall  pressure 
coefficient  along  the  circumferential  direction  of  a  cylinder  at  Re  =  40  is  shown 
in  Figure  5.1.  On  the  average,  the  pressure  field  was  correct,  but  the  oscillations 
were  very  severe.  The  oscillations  were  a  consequence  of  the  use  of  a  non-staggered 
grid.  At  this  point,  the  present  pressure  correction  formulation  was  chosen.  Using 
a  similar  pressure  correction  formulation,  Rider  [108]  found  omitting  the  pressure 
term  in  the  auxiliary  velocity  equation  produces  a  formulation  which  has  better 
error  characteristics,  without  a  loss  of  accuracy,  compared  to  a  pressure  correction 
method  which  accumulates  the  pressure  corrections. 
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5.2  Diagonal  Dominance 


The  form  of  the  matrix  corresponding  to  the  resulting  discrete  system  following 
discretization  of  the  Poisson  equation  is  critical  to  convergence  of  iterative  methods 
employed  to  invert  it.  Diagonal  dominance  is  the  key  property  of  the  matrix  related 
to  convergence  of  an  iterative  method.  Matrix  A  of  order  n  is  (strictly  row)  diagonal 
dominant  [35]  if 


\aii\  ^  ^  v  \aij  1?  ^  ~~  (5.10) 

Diagonal  dominance  of  the  matrix  is  dependent  upon  the  procedure  used  to 
generate  the  matrix.  To  examine  this  property,  we  used  a  finite  difference,  finite 
element,  and  finite  volume  technique  to  generate  the  Poisson  matrix.  To  enable 
the  use  of  the  finite  difference  technique,  a  structured  prismatic  grid  was  used  in 
the  analysis,  similar  to  the  grid  shown  in  Figure  6.5.  The  Gauss-Seidel  method 
was  used  to  solve  the  system  of  equations  for  all  three  cases.  The  Poisson  equation 
convergence  criteria  was  also  identical  for  all  three  cases.  The  quasi- two-dimensional 
flow  around  a  cylinder  at  a  Reynolds  number  of  40  was  the  chosen  test  case. 

An  indication  of  the  level  of  diagonal  dominance  of  the  matrix  can  be  seen 
by  examining  the  ratio 


Ri  = 


jo»  I 

\aij\ 


(5.11) 


If  each  row  i  has  a  diagonal  dominance  ratio  greater  than  1.0,  the  matrix  is  diagonally 
dominant.  The  diagonal  dominance  ratio  of  the  resulting  matrix  for  a  specific  row 
is  examined  in  Table  5.1  for  the  finite  difference,  finite  element,  and  the  finite 
volume  technique.  The  finite  difference  method  matrix  had  a  ratio  of  1.0,  which  is 
(marginally)  diagonally  dominant.  The  finite  element  method  produced  a  matrix 
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Table  5.1:  Indicators  of  the  level  of  diagonal  dominance  for  Poisson  matrix  formu¬ 
lations.  A  ratio  greater  than  1.0  signifies  diagonal  dominance. 


Finite  Difference 

Finite  Element 

Finite  Volume 

Diagonal 
Dominance 
Ratio,  R{ 

1.0 

0.91 

0.22 

with  a  ratio  of  0.91,  a  slightly  non-diagonally  dominant  matrix.  The  finite  volume 
discretization  produced  a  matrix  with  a  ratio  of  0.22,  a  strongly  non-diagonally 
dominant  matrix. 

A  comparison  of  the  number  of  Gauss-Seidel  iterations  required  for  the  con¬ 
vergence  for  the  three  cases  is  shown  in  Figure  5.2.  For  the  finite  difference  (diag¬ 
onally  dominant)  matrix,  the  number  of  Gauss-Seidel  iterations  at  each  time  step 
was  initially  high,  but  quickly  decreased  to  one  Gauss-Seidel  iteration  per  time  step. 
The  finite  element  (slightly  non-diagonally  dominant)  matrix  had  a  similar  conver¬ 
gence  behavior.  For  the  finite  volume  (strongly  non-diagonally  dominant)  matrix, 
convergence  was  not  obtained.  The  momentum  residuals  were  also  compared  and 
are  shown  in  Figure  5.3.  The  solution  using  the  finite  volume  method  (strongly 
non-diagonally  dominant)  never  converged.  The  convergence  behavior  of  the  finite 
difference  and  finite  element  methods  are  very  similar. 

For  semi-unstructured  and  unstructured  meshes,  the  finite  difference  tech¬ 
nique  is  unsuitable  because  it  cannot  handle  the  arbitrariness  of  node  placement. 
The  loss  of  diagonal  dominance  makes  the  standard  finite  volume  technique  for  the 
Poisson  equation  inadequate.  For  this  reason,  a  finite  element  technique  was  used 
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Figure  5.2:  Comparison  of  number  of  Gauss-Seidel  iterations  needed  for  convergence 
of  pressure  Poisson  equation. 

- Finite  difference  (diagonally  dominant)  matrix 

—  Finite  element  (slightly  non- diagonally  dominant)  matrix 
- Finite  volume  (strongly  non- diagonally  dominant)  matrix 
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Figure  5.3:  Comparison  of  maximum  momentum  residual  behavior. 

- Finite  difference  (diagonally  dominant)  matrix 

—  Finite  element  (slightly  non-diagonally  dominant)  matrix 
- Finite  volume  (strongly  non-diagonally  dominant)  matrix 
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for  solution  of  the  Poisson  equation. 

This  difficulty  with  the  formation  of  non-diagonally  dominant  matrices  is 
usually  not  an  issue  for  staggered  grids  used  with  the  finite  volume  technique.  This  is 
related  to  the  number  of  neighboring  nodes  needed  to  compute  the  Poisson  equation 
for  each  node.  When  the  finite  volume  technique  is  applied  to  a  node  in  a  non- 
staggered  hexahedral  mesh,  the  pressure  values  at  all  the  neighboring  nodes  are 
required.  The  resulting  Poisson  matrix  will  have  26  off-diagonal  entries  in  each 
row  of  the  matrix,  and  is  usually  not  row  diagonally  dominant.  When  the  finite 
volume  technique  is  applied  to  node  in  a  pressure-staggered  (velocity  is  defined  at 
the  nodes,  pressure  is  defined  at  the  cell  centers)  hexahedral  mesh,  the  pressure 
values  at  the  neighboring  cell  centers  is  required.  The  resulting  Poisson  matrix  will 
have  only  8  off-diagonal  entries  in  each  row  of  the  matrix.  The  matrix  resulting 
from  a  hexahedral  staggered  mesh  is  diagonally  dominant  [70]. 


60 


5.3  Boundary  Conditions  for  the  Poisson  Equation 


The  pressure  correction  equation  is  a  Poisson  equation  and  can  be  written  in  the 
following  form. 

A  <j>  =  -/  (5.12) 

The  value  of  the  function  <f>  on  a  boundary,  can  be  specified  in  two  different  ways. 
The  first  type  of  boundary  condition  specifies  the  value  of  the  function  <f>  on  the 
boundary  surface  (Dirichlet  boundary  conditions): 

<t>  =  g{s)  (5.13) 

The  second  type  of  boundary  condition  specifies  the  value  of  the  normal  derivative 
of  <f>  on  the  boundary  surface  (Neumann  boundary  conditions): 

tn  =  S{S)  (5'14) 

Dirichlet  boundary  conditions  with  the  Poisson  equation  have  been  used 
sparingly  in  the  literature.  Improper  definition  of  the  boundary  condition  may  create 
local  mass  sources  and  create  sinks  to  appear  [45].  Most  of  the  methods  published 
in  the  literature  use  Neumann  boundary  conditions  for  the  Poisson  equation.  The 
choice  of  using  Dirichlet  or  Neumann  boundary  conditions  has  been  discussed  in  the 
literature  [1,  3,  43]. 

In  the  present  work,  severe  pressure  oscillations  resulted  when  <f>  was  speci¬ 
fied.  Moving  the  far  field  boundary  further  away  from  the  body  did  not  result  in  a 
decrease  in  the  pressure  oscillations.  Since  severe  oscillations  are  often  an  indication 
of  unsuitable  boundary  conditions  [44],  Neumann  boundary  conditions  were  used. 

For  the  Neumann  problem,  g  may  not  be  chosen  independently  of  /,  but 
rather  as  a  consequence  of  Green’s  first  identity.  The  identity  states  /  and  g  must 
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satisfy  a  compatibility  condition, 

f  fdil  =  [  g(s)ds  (5.15) 

Ju  JdU 

if  a  solution  is  to  exist  [3,  22,  42].  A  property  of  the  Neumann  problem  is  the 
non-uniqueness  of  the  solution.  Although  the  shape  of  the  solution  is  unique,  the 
location  of  this  solution  is  not  fixed  in  space,  unlike  the  Dirichlet  problem.  There 
are  an  infinite  number  of  solutions  each  differing  from  the  others  by  an  arbitrary 
additive  constant. 

One  method  developed  to  eliminate  the  floating  of  the  solution  is  to  choose 
one  point  in  the  flow  field  and  assign  to  it  a  specific  value,  essentially  forcing  one 
point  to  have  Dirichlet  boundary  conditions  [134,  135].  The  choice  of  selecting  an 
interior  point  has  also  been  studied  [56].  However,  the  selection  of  a  Dirichlet  point 
may  cause  the  formation  of  a  singularity  at  that  point.  Failure  to  recompute  the 
value  at  this  point  results  in  the  partial  differential  equation  not  being  satisfied  at 
this  point.  Thus  the  link  with  any  solution  is  lost  at  this  point  [38].  Accumulation  of 
nearly  all  of  the  round-off  error  may  also  occur  at  this  point,  sending  perturbations 
throughout  the  entire  flow  field  [56]. 

A  second  method  developed  to  eliminate  the  floating  of  the  solution  also 
selects  one  point  in  the  flow  field.  However,  the  solution  is  computed  for  all  points, 
even  the  selected  point.  The  resultant  approximation  to  the  solution  is  shifted  by 
subtracting  the  value  of  the  solution  at  the  selected  point  from  all  the  points  in  the 
flow  field  [38].  The  use  of  this  method  results  in  successive  iterations  approximating 
the  same  solution. 

A  third  method  developed  to  eliminate  the  floating  of  the  solution  uses  the 
compatibility  condition,  equation  (5.15),  as  an  integral  constraint  on  the  pressure 
[22,  69,  103].  This  method  produces  acceptable  results,  but  is  computationally  more 
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expensive  than  the  second  method. 

In  the  present  work,  the  second  method  is  used  to  eliminate  the  floating  of 
the  <f>  values,  and  subsequently  the  pressure  values.  One  point  in  the  far  field  is 
chosen,  and  the  resulting  pressure  value  is  subtracted  from  the  entire  pressure  field 
each  time  the  pressure  is  corrected. 
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Figure  5.4:  Placement  of  primitive  variables  on  a  non-st aggered  grid. 


5.4  Non-staggered  Grid  versus  Staggered  Grid 


On  a  non-staggered  grid,  values  of  all  the  primitive  variables  are  defined  at  each 
node.  Figure  5.4  illustrates  a  non-staggered  grid  in  two-dimensions.  Figure  5.5 
illustrates  a  one-dimensional  uniform  domain  used  to  calculate  the  pressure  gradient 
at  node  i  using  the  finite  volume  approach  on  a  non-staggered  grid.  The  pressure 
gradient  at  node  i  is  found  from  the  difference  of  the  pressure  value  at  faces  e  and  w. 
Since  the  pressure  is  defined  at  the  nodes  only  (non-staggered  grid)  the  value  of  the 
pressure  at  each  face  is  found  by  linearly  interpolating  the  values  of  the  pressure  at 
the  two  adjacent  nodes  as  shown  in  equation  (5.16).  However,  the  equation  simplifies 
and  the  resulting  pressure  gradient  term  will  contain  the  pressure  difference  between 
two  alternate  grid  points,  and  not  between  adjacent  ones. 


m  =  Pe~Pw 
dx  Ax 


Ax 


Pj+ 1  —  Pj-i 

2Ax 


(5.16) 
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Figure  5.5:  One-dimensional  domain  for  pressure  gradient  calculation  on  a  non- 
staggered  grid. 

Comparison  of  the  smooth  and  odd-even  pressure  fields  on  the  one-  dimen¬ 
sional  domain  are  shown  in  Figures  5.6  and  5.7.  The  resulting  pressure  gradient 
calculation  from  equation  (5.16)  produces  the  same  result  for  both  pressure  fields. 
A  consequence  of  central  differencing  on  a  non-staggered  grid  is  that  a  wavy  pres¬ 
sure  field  will  be  “felt”  like  a  uniform  pressure  field  by  the  momentum  equation.  If 
oscillatory  pressure  fields  arise  during  the  iterative  pressure  correction  procedure, 
the  oscillations  would  be  preserved  until  convergence  since  the  momentum  equations 
would  be  oblivious  to  their  presence. 

This  difficulty  can  be  eliminated  with  the  use  of  staggered  grids  [100].  On  a 
staggered  grid,  the  values  of  all  the  primitive  variables  are  defined  at  different  nodes. 
Figure  5.8  illustrates  one  example  of  a  staggered  grid.  The  staggering  concept  may 
be  different  than  shown  in  Figure  5.8.  The  pressure  staggered  grid  [70]  defines  the 
velocities  at  cell  nodes  and  the  pressure  is  defined  at  the  cell  center.  The  arbitrary 
Lagrangian-Eulerian  (ALE)  procedure  [27,  54,  55]  defines  the  pressure  at  the  nodes 
and  velocity  components  at  the  cell  center. 

In  the  past,  most  incompressible  flow  solvers  used  staggered  grids  [47,  48, 
52,  67,  70,  72,  123].  A  one- dimensional  staggered  grid  analogy  to  Figure  5.5  is 
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Figure  5.6:  Uniform  pressure  field  in  a  one-dimensional  domain. 


Figure  5.7:  Odd-even  pressure  field  in  a  one-dimensional  domain. 
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Figure  5.8:  Placement  of  primitive  variables  on  a  staggered  grid. 
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Figure  5.9:  One-dimensional  domain  for  pressure  gradient  calculation  on  a  staggered 
grid. 


illustrated  in  Figure  5.9.  The  figure  illustrates  a  one-dimensional  uniform  domain 
used  to  calculate  the  pressure  gradient  at  node  i  using  the  finite  volume  approach  on 
a  staggered  grid.  The  advantage  of  the  staggered  grid  is  that  the  pressure  gradient 
at  node  i  is  calculated  from  the  pressure  at  the  faces  e  and  w,  which  are  known. 
Therefore,  the  pressure  gradient  at  node  i  is  calculated  from  the  two  adjacent  grid 
points,  and  will  not  permit  oscillatory  pressure  fields.  A  disadvantage  of  the  use  of 
staggered  grids  is  that  each  staggered  variable  requires  its  own  control  volume.  For 
the  staggered  grid  shown  in  Figure  5.8,  three  control  volumes  would  be  required. 
Two  control  volumes  for  the  u  and  v  velocities  and  one  control  volume  for  the 
pressure  would  necessitate  an  increase  in  memory  storage. 

The  non-staggered  grid  requires  only  one  control  volume  for  all  the  vari¬ 
ables.  This  is  an  advantage  in  terms  of  memory  storage  and  coding  complex¬ 
ity.  As  the  popularity  of  unstructured  grids  rises,  the  use  of  non-staggered  grids 
[10,  49,  106,  122,  126,  130]  has  become  more  prevalent.  The  present  work  uses  hy¬ 
brid  (prismatic/ tetrahedral)  meshes  which  led  to  the  choice  of  non-staggered  grid. 
Care  must  be  taken  to  ensure  an  oscillatory  pressure  field  does  not  develop  during 
the  iterative  pressure  correction  process. 

A  comparison  of  solutions  and  convergence  properties  of  a  finite  volume  nu- 
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merical  method  with  staggered  and  non-staggered  grids  [102]  concluded  the  non- 
staggered  version  had  no  disadvantages  in  performance  relative  to  the  staggered  ver¬ 
sion.  In  fact,  treatment  of  boundary  conditions  and  the  implementation  of  higher 
order  differencing  schemes  was  simpler  for  the  non-staggered  version.  In  some  cases, 
the  non-staggered  version  displayed  faster  convergence.  Convergence  acceleration 
techniques,  like  the  multigrid  technique,  are  also  easier  to  apply  to  the  non-staggered 
grid  arrangement. 

A  disadvantage  to  non-staggered  methods  not  mentioned  by  [102],  but  found 
in  all  non-staggered  methods  [10,  102,  106,  124]  is  a  non-zero  value  of  the  divergence 
of  the  velocity.  Further  explanation  of  this  phenomenon  is  included  in  the  subsection 
on  Discrete  Divergence. 
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5.5  Discrete  Divergence 


A  potentially  serious  consequence  to  numerically  modeling  a  flow  field  is  errors  in¬ 
troduced  to  the  solution  due  to  the  discretization  of  the  governing  equations.  For 
the  present  work,  the  momentum  and  pressure  Poisson  equations  are  the  governing 
equations  of  interest.  These  equations  can  be  discretized  independently  of  one  an¬ 
other  or  the  discretization  of  one  equation  can  be  developed  using  the  other  equation 
in  a  consistent  formulation.  Using  a  methodology  similar  to  [124],  the  discretization 
procedure  of  the  present  pressure  correction  formulation  for  one  dimension  in  finite 
difference  form  is  explained.  The  1-D  momentum  equation  is  considered. 


du  _ 

dt  dx 


(5.17) 


where 

1  d2u 
dx  Re  dx 2 

The  1-D  auxiliary  velocity  vector  is 


(5.18) 


(5.19) 


5.5.1  Consistency  of  Discretized  Poisson  Equation 

The  discrete  Poisson  equation  will  be  derived  using  two  methods,  to  determine  the 
consistency  of  the  formulation.  The  methods  are: 

1.  The  Poisson  equation  is  derived  using  the  momentum  and  continuity  equa¬ 
tions.  The  final  form  of  the  Poisson  equation  is  then  discretized. 

2.  The  discrete  Poisson  equation  is  derived  used  the  discrete  momentum  and  the 
discrete  continuity  equations. 
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The  resulting  discrete  Poisson  equations  will  be  compared  for  consistency  using 
both  staggered  and  non-staggered  grids.  Central  differencing  will  be  used  in  the 
discretizations. 


Staggered  Grid 

In  the  first  method,  the  Poisson  equation  is  derived  using  the  momentum  and  con¬ 
tinuity  equations.  The  Poisson  equation  in  one  dimension,  which  was  derived  pre¬ 
viously  and  originally  defined  in  equation  (2.36)  is  repeated  here  for  clarity 


A  cf>  =  V  •  uf 


(5.20) 


and  rewritten  in  terms  of  the  pressure  variable. 

V-i? 


A  P  = 


At 


(5.21) 


Using  a  staggered  grid  where  the  pressure  is  defined  at  the  nodes  (i-l,i,i+l)  and 
the  velocity  is  staggered  between  the  nodes  (i+1/2  and  i-1/2),  the  discrete  auxiliary 
vector  at  locations  i+1/2  and  i-1/2  become: 


Ui+l/2  —  <+1/2  ~  1/2 


(5.22) 


(5.23) 


<-t/ 2  =  <-1/2  -  AtfT-i/a 

Substitution  of  equations  (5.22)  and  (5.23)  into  equation  (5.21),  produces  the  dis¬ 
crete  Poisson  equation. 

(Fg?-2Pp'  +  rtf\  +  n\  _  i_  Knn-'Un 


(5.24) 


y  (Ax)2  J  '  ^  Ax  J  A t  y  Ax 

Applying  the  second  method  on  a  staggered  grid,  the  discrete  Poisson  equa¬ 
tion  is  derived  using  the  discrete  continuity  and  discrete  momentum  equations.  The 
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discrete  momentum  equation  at  locations  i+1/2  and  i-1/2  are: 

-  At£”+1/2 


<+1/2  =  <+1/2  “  At 


f  pn+\  __  pn+1  ^ 

A# 


<-1/2  =  <-1/2  "  At 


'  jpn+1  __  j^n+l ' 


Ax 


At^JLx/2 


The  discretized  continuity  equation  is: 


7/"+1  -  ?/n+1 

»+l/2  »— 1/2  _ 


Ax 


=  0 


(5.25) 

(5.26) 

(5.27) 


Substituting  the  discrete  momentum  equations  (5.25)  and  (5.26)  into  the  discrete 
continuity  equation  (5.27),  it  is  obtained: 


(p?A'-u?+'  +  i?+')  ((U/2-dw)  =  j_ 

Ax 


(5.28) 


\  (Ax)2  J  *  \  A x  J  At 

Comparisons  of  equations  (5.24)  and  (5.28),  show  the  resulting  discrete  Poisson 
equation  obtained  by  both  methods  are  identical  on  a  staggered  grid.  The  staggered 
grid  formulation  is  consistent. 


Non-staggered  Grid 

In  a  non-staggered  grid,  the  velocity  and  pressure  are  defined  at  the  nodes  (i-l,i,i+l). 
Using  the  first  method  to  calculate  the  discrete  Poisson  equation  requires  the  discrete 
auxiliary  vector  at  nodes  i+1  and  i-1: 


=  <+i  -  A^r+i 

uT.Y  =  <_x  -  Atff_i 


(5.29) 

(5.30) 


Substitution  of  equations  (5.29)  and  (5.30)  into  the  Poisson  equation  (5.21)  produces 
the  discrete  Poisson  equation. 
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Applying  the  second  method  on  a  non-staggered  grid,  the  discrete  Poisson 
equation  is  derived  using  the  discrete  continuity  and  discrete  momentum  equations. 
The  discrete  momentum  equations  at  nodes  i-fl  and  i-1  are: 


(pn+1  pn +1  \ 

2/S.x  ) 


(  pn4-l  _  pn4-l  ’ 

tlft1  =  «?_!  -At  *  *-2 


2Az 


-  Ater+i 


The  discretized  continuity  equation  is: 


1/"+1  _ 

Mi+i  ui-i 

2Ax 


=  0 


(5.32) 

(5.33) 


(5.34) 


Substituting  the  discrete  momentum  equations  (5.32)  and  (5.33)  into  the  discrete 
continuity  (5.34),  it  is  obtained: 


Comparisons  of  equations  (5.31)  and  (5.35),  show  the  resulting  discrete  Poisson 
equations  obtained  by  the  two  methods  are  not  identical  on  a  non-staggered  grid. 
The  non-staggered  grid  formulation  is  not  consistent. 


5.5.2  Error  of  Discrete  Continuity 

The  discrete  Poisson  equation,  (5.28),  using  staggered  grids  is  obtained  using  a 
consistent  formulation.  It  is  observed  that  the  right  hand  side  of  the  equation  is  the 
discrete  continuity  equation,  while  the  left  hand  side  is  the  difference  of  the  residual 
of  the  momentum  equation  at  two  different  nodes.  At  steady  state,  the  residual  of 
the  momentum  equation  at  all  nodes  is  zero.  If  the  left  hand  side  of  the  equation 
is  driven  to  zero,  so  will  the  right  hand  side.  Therefore,  for  a  staggered  grid,  the 
discrete  continuity  equation  will  approach  machine  zero. 
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The  discrete  Poisson  equations,  (5.31)  and  (5.35),  using  non-staggered  grids 
are  not  consistent.  The  difference  between  equations  (5.31)  and  (5.35)  is  the  error  in 
the  discrete  continuity.  To  examine  this  error,  the  pressure  term  in  equation  (5.35) 


will  be  subtracted  from  and  added  to  the  left  hand  side  of  equation  (5.31). 

(P?+i  -  2 P[^  +  P-X  \  (P$  -  2 P[+1  +  P?_+l 


(Ax)2 


W 


(2Ax)2 


+ 


( Pi+2  -  2Pll+l  +  P?-2  \  .  ftf+1  ~  tf-1  \  1  (<+l  -  <-!  \  ,E  ^ 

{ - (2A^ - J  +  {  2Ax  )  =  A 1  {  2Ax  )  (5’36) 


The  third  and  fourth  terms  on  the  left  hand  side  of  the  equation  are  the  difference  of 
the  residuals  of  the  momentum  equation  at  two  different  nodes  on  a  non-staggered 
grid.  At  steady  state,  those  residuals  will  be  zero.  Removing  the  third  and  fourth 
terms  from  the  equation  produces 

'  Pi+2  ~  2-P"+1  +  Pi-2 

(2Ax)2  j  A f  V  2Ax  ) 

(5.37) 


(P&1  -2P?+1  +P?_+1\  L 

\  (A^)2  J  {' 


when  rearranged  becomes 


PI a1  +  -  6ff(1  +  4ig;i  -  jy_V \  ... 

4(Ax)2  )  ~  (  2Ax  J  (5'38) 


which  can  be  rewritten  in  a  more  concise  form  as 


At  d4P 
A(Ax)2  dx 4 


(5.39) 


The  right  hand  side  of  this  equation  is  the  discrete  continuity.  The  left  hand 
side  of  the  equation  is  the  error  in  the  discrete  continuity  with  non-staggered  grids 
at  steady  state.  The  discrete  continuity  residual  approaches  a  term  proportional 
to  the  fourth  derivative  of  the  pressure,  the  time  step  and  the  cube  of  the  grid 
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spacing,  which  is  the  same  error  found  by  [124]  using  a  different  pressure  correction 
formulation  on  a  non-staggered  grid.  Since  the  time  step  size,  A/,  is  proportional  to 
the  spatial  step  size,  Ax,  the  order  of  the  discrete  divergence  error  is  (Ax)3,  which 
is  smaller  than  the  spatial  discretization  error  (Ax)2. 

For  non-staggered  grids,  the  Poisson  equation  is  discretized  inconsistently 
relative  to  the  momentum  equations.  This  results  in  the  discrete  divergence  of  the 
velocity  not  being  driven  to  machine  zero.  It  has  been  shown  by  others  the  discrete 
divergence  can  not  be  driven  to  machine  zero,  and  at  the  same  time,  still  maintain 
a  smooth  pressure  field  [124,  125,  127].  To  obtain  a  smooth  pressure  field  on  a 
non-staggered  grid,  [125]  adds  an  artificial  mass  source  term  to  the  right  hand  side 
of  the  discrete  continuity  equation.  To  enforce  mass  conservation  for  internal  flows, 
[17,  122]  imposes  the  known  mass  flow  rate  at  the  inlet  and  outlet  boundaries. 

Several  papers  in  the  literature  which  use  similar  methods  (projection  meth¬ 
ods  on  non-staggered  grids)  have  claimed  their  methods  conserve  mass  exactly 
[107,  147].  These  statements  have  been  described  as  misleading  by  [108].  Their 
cell  centered  methods  interpolate  the  velocities  to  the  cell  edges  and  the  divergence 
is  constructed  from  these  values.  The  discrete  divergence  for  the  edge-centered  di¬ 
vergences  will  exactly  satisfy  continuity.  However,  the  velocities  at  the  cell  centers 
will  not  in  general  satisfy  continuity  exactly.  Because  the  cell  centered  data  are 
advanced  in  time,  [108]  feels  it  is  invalid  to  state  those  methods  conserve  mass 
exactly.  The  momentum  interpolation  method  also  has  difficulty  with  the  velocities 
not  satisfying  the  continuity  equation  [54]. 
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5.6  Analytic  Field  Function  Testing 


Critical  subroutines  of  the  solver  are  tested  using  analytic  field  functions.  The  di-  > 
vergence  calculation  is  tested  using  a  known  analytic  function.  The  finite  volume 
divergence  calculation  is  very  similar  to  the  calculation  used  in  the  convective  and 
gradient  calculations,  so  the  testing  of  these  subroutines  will  not  be  undertaken. 
Since  the  function  being  operated  on  is  a  known  analytic  function,  the  resulting 
solution  is  also  a  known  analytic  function.  This  analytic  solution  is  then  compared 
with  the  solution  obtained  from  the  numerical  calculation.  This  testing  of  the  solver 
using  analytic  field  functions  not  only  shows  the  basic  accuracy  of  the  solver  sub¬ 
routines,  but  also  gives  a  quantitative  evaluation  of  the  suitability  of  the  grid.  The 
location  of  high  error,  during  the  analytic  test,  can  be  identified,  but  if  the  errors 
occur  in  locations  of  the  grid  where  changes  in  the  flow  variables  are  insignificant, 
the  errors  would  not  occur  during  an  actual  flow  simulation  and  the  grid  is  suitable 
for  a  simulation.  The  percent  error  of  the  numerical  solution  is  calculated  using 

ny  t-,  ,  ( numerical  —  analytic\  H  ^  . 

%Error  =  abs  - - — : - - —  *  100  (5.40) 

\  analytic  J 

In  addition,  the  effect  of  grid  size  on  the  accuracy  of  divergence  and  viscous  calcu¬ 
lations  is  examined  with  three  successively  finer  grids. 

5.6.1  Analytic  Field  Test  of  Sphere  Mesh 

An  all  prismatic  mesh  and  a  hybrid  mesh  for  a  sphere  are  used  to  demonstrate  the 
accuracy  of  the  numerical  calculations.  Grid  parameters  of  the  meshes  are  shown  in 
Table  5.2. 
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Table  5.2:  Grid  parameters  used  for  sphere  meshes  for  analytic  field  function  anal¬ 
ysis. 


Type 

of 

Grid 

Number  of 
Surface 
Faces 

Number  of 
Prism 
Cells 

Number  of 
Tetrahedral 
Cells 

Total 
Number 
of  Nodes 

Radial 

Step 

Size 

Prism 

1346 

87490 

0 

46530 

0.01 

Hybrid 

1346 

36342 

160054 

49195 

0.01 

Divergence  of  a  Vector 

The  divergence  calculation  is  tested  using  a  vector  which  is  a  linear  function  in  x,y,z. 
The  analytic  linear  vector  is 

U  =  (x  +  y  +  z)i  +  (x  +  y  +  z)j  +  (x  +  y  +  z)k  (5-41) 

The  analytic  solution  is 

V  •  U  =  3.0  (5.42) 

The  numerical  errors  of  the  divergence  calculation  of  a  linear  function,  using 
the  prismatic  mesh  and  the  hybrid  mesh,  are  shown  in  Figure  5.10.  The  percentage 
of  interior  nodes  is  plotted  against  the  percent  error  in  increments  of  0.1  percent. 
It  is  observed  that  the  majority  of  the  interior  nodes  in  the  prismatic  mesh  have  an 
error  of  less  than  1.0%.  Maximum  error  for  this  mesh  was  1.8%.  The  majority  of 
the  nodes  in  the  hybrid  mesh  have  errors  less  than  0.1%.  Nodes  contributing  the 
highest  errors  are  located  at  the  interface  of  the  prismatic  and  tetrahedral  regions. 


76 


%  Error 

Figure  5.10:  Percent  error  versus  percent  interior  nodes  for  the  divergence  calcula¬ 
tion  of  a  linear  analytic  function  around  a  sphere. 

- Prismatic  mesh. 

- Hybrid  mesh. 

5.6.2  Order  of  Accuracy  Evaluation 

The  effect  of  grid  size  on  accuracy  is  demonstrated  with  three  semi-unstructured 
prismatic  meshes  around  a  sphere.  Only  one  quadrant  of  the  sphere  is  modeled. 
The  medium  mesh  is  created  by  adding  a  radial  node  between  every  radial  node 
in  the  coarse  mesh  and  adding  a  node  in  the  center  of  the  triangular  faces  in  the 
prismatic  cells.  One  coarse  prism  cell  is  subdivided  into  eight  prism  cells  in  order 
to  construct  the  medium  mesh.  The  fine  mesh  is  created  from  the  medium  mesh 
using  the  same  refinement  principles.  Specific  parameters  of  the  meshes  are  shown 
in  Table  5.3.  The  root  mean  square  (rms)  error  is  defined  as 

RMS  Error  =  ^numerical  -  anatyUorf 

total  interior  nodes 

The  analytic  function  used  in  the  divergence  calculation  was  shown  in  equa- 
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Table  5.3:  Grid  parameters  of  sphere  meshes  for  analytic  field  function  grid  fineness 
study. 


Type 

of 

Grid 

Number  of 
Radial 
Nodes 

Number  of 
Surface 
Nodes 

Number  of 
Surface 
Faces 

Total 
Number 
of  Nodes 

Radial 

Step 

Size 

Coarse 

51 

87 

139 

4437 

0.08 

Medium 

101 

311 

556 

31411 

0.04 

Fine 

201 

1176 

2224 

236376 

0.02 

tion  (5.41),  while  the  analytic  function  used  for  the  viscous  calculation  is: 

u  —  x2i  +  y2j  +  z2k  (5.44) 

The  analytic  solution  is 

V2u  =  6.0  (5.45) 

The  RMS  error  of  the  divergence  and  viscous  calculations  using  the  coarse, 
medium,  and  fine  grids  are  shown  in  Figure  5.11.  The  slope  of  the  log-log  plot  shows 
both  calculations  are  more  accurate  as  the  mesh  is  refined  and  have  second  order 
accuracy. 
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Chapter  6 


Validation  and  Application  of 
the  Adaptive  Grid  Solver 


Validation  of  the  three-dimensional  hybrid  Navier-Stokes  solver  for  incompressible 
flows  is  demonstrated  in  a  piecewise  fashion.  First,  stability  of  the  solver  at  a  high 
Reynolds  number  is  examined.  Spatial  accuracy  of  the  solver  is  investigated  on 
steady  quasi-two-dimensional  cylinder  flow  using  two  different  prismatic  meshes. 
Time  accuracy  is  examined  with  an  impulsive  start  of  flow  over  a  cylinder  using  a 
prismatic  mesh.  Three-dimensional  aspects  of  the  solver  are  validated  using  driven 
cubic  cavity  flow,  as  well  as  flow  over  a  sphere  using  prismatic  meshes  at  several 
Reynolds  numbers.  The  effectiveness  of  prismatic  adaptation  is  demonstrated  with 
cubic  cavity  flow.  Equation  and  hybrid  grid  adaptation  are  also  illustrated  with 
steady  flow  over  a  sphere.  The  effectiveness  of  hybrid  adaptation  is  also  demon¬ 
strated  with  flow  over  tandem  spheres. 
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6.1  Flow  over  a  Flat  Plate  at  a  High  Reynolds  Number 


To  examine  the  stability  of  the  algorithm  at  high  Reynolds  numbers,  the  flow  over 
a  flat  plate  is  computed  at  Re  =  106,  based  on  the  length  of  the  plate.  Grid 
parameters  are  shown  in  Table  6.1  while  boundary  conditions  are  shown  in  Figure 
6.1.  The  computational  domain  begins  at  x=1.0  with  inlet  conditions  of  a  Blasius 
profile.  The  artificial  dissipation  factor,  (74,  was  equal  to  10”6  and  the  CFL  factor 
was  0.01.  The  signature  of  the  semi-unstructured  prismatic  mesh  on  the  plate,  as 
well  as  on  the  other  boundaries,  is  shown  in  Figure  6.2.  The  plate  is  tessellated  with 
triangles,  while  the  side  boundaries  are  tessellated  with  the  quadrilateral  faces  of 
the  prisms. 

The  u  velocity  profile  at  the  center  of  the  plate  is  shown  in  Figure  6.3.  The 
numerical  results  coincide  with  the  Blasius  velocity  profile  [120].  The  convergence 
history  of  the  simulation  in  terms  of  the  maximum  residual  in  the  x-momentum 
equation  is  shown  in  Figure  6.4.  The  maximum  residual  is  defined  as  the  maximum 
of  | un+1  —  un |  over  all  nodes.  From  these  results,  the  following  conclusions  can  be 
drawn: 

1.  The  pressure  correction  scheme  with  a  non-staggered  grid  converges  well  and 
is  stable  for  flow  over  a  flat  plate  at  a  high  Reynolds  number. 

2.  The  incompressible  three-dimensional  solver,  using  a  semi-unstructured  pris¬ 
matic  mesh,  yields  an  accurate  solution  when  compared  to  analytical  results. 
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Table  6.1:  Grid  parameters  used  for  flat  plate  flow  at  Re=  106. 


Number  of 
Faces  on  Wall 

Grid  Spacing 
at  Wall 

Grid  Stretching 
Factor 

Number  of 
Prism  Layers 

Number  of 
Nodes 

1774 

1.125 

70 

67308 

Far  Field 

du/dz  =  v  =  dw/dz  =  0 
d<))  /dz  =  0 


Outlet 


du/dx  =  0 
dw/dx  =  0 
v  =  0 
d())/dx  =  0 


u  =  v  =  w  =  0 

d())/dz  =  0  x-z  Symmetry  Planes 

du/dy  =  v  =  dw/dy  =  0 
d<j)  /dy  =  0 


Figure  6.1:  Boundary  conditions  for  flat  plate  simulation  at  Re  =  106. 
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Figure  6.2:  Prismatic  mesh  with  unstructured  surface  tessellation  for  flat  plate  at 
Re  =  106. 


0.020 


0.015 


Figure  6.3:  U  velocity  profile  at  the  center  of  a  flat  plate  (Re  =  106). 
o  Blasius  result  —  Present  work 


Figure  6.4:  Convergence  of  the  maximum  residual  in  the  x-momentum  equation  for 
a  flat  plate  (Re  =  106). 
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6.2  Steady  Flow  Around  a  Cylinder  using  a  Prismatic 


Mesh 

Steady  quasi-two-dimensional  flow  around  a  circular  cylinder  at  a  Reynolds  number 
of  40,  based  on  the  cylinder  diameter,  is  used  to  examine  the  spatial  accuracy  of 
the  solver  on  a  prismatic  mesh.  At  this  Reynolds  number,  the  flow  is  attached  to 
the  cylinder  over  most  of  its  surface,  with  an  attached  separation  bubble  behind  the 
cylinder. 

Two  meshes  were  used  to  demonstrate  the  validity  of  the  solver.  The  first 
mesh  has  a  “structured”  surface  triangulation,  while  the  second  mesh  has  an  “un¬ 
structured”  surface  triangulation.  The  first  mesh  is  shown  in  Figure  6.5  and  is 
generated  by  subdivision  of  quadrilateral  faces  along  one  of  their  diagonals.  The 
second  “unstructured”  mesh  is  shown  in  Figure  6.6  which  is  generated  via  an  ad¬ 
vancing  front-type  of  surface  grid  generation.  The  cylinder  axial  length  is  three 
times  the  diameter.  Grids  of  varying  fineness  were  tested  until  the  results  were  no 
longer  mesh  dependent.  The  grid  parameters  used  for  this  case  are  shown  in  Table 
6.2.  Symmetry  planes  were  defined  at  both  ends  of  the  cylinder.  Location  of  the 
far  field  boundary  is  at  approximately  21  diameters.  Surface  triangular  faces  were 
stretched  in  the  axial  direction  since  the  flow  does  not  vary  along  the  axis  of  the 
quasi  two-dimensional  cylinder.  The  artificial  dissipation  factor,  <j4,  was  equal  to 
10“5  and  the  CFL  factor  was  0.1.  Results  from  the  three-dimensional  solver  show 
the  flow  is  basically  two-dimensional. 

Particle  paths  for  one  cross-sectional  plane  of  the  cylinder  are  shown  in  Figure 
6.7.  These  paths  were  obtained  from  the  solution  obtained  from  the  “unstructured” 
mesh.  The  attached  flow  and  separation  bubble  are  clearly  shown. 
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Figure  6.5:  Partial  prismatic  mesh  for  circular  cylinder  at  Re  =  40  with  “structured” 
surface  triangulation. 


Figure  6.6:  Partial  prismatic  mesh  for  circular  cylinder  at  Re  =  40  with  “unstruc¬ 
tured”  surface  triangulation. 
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Table  6.2:  Grid  parameters  used  for  circular  cylinder  flow  at  Re=40. 


Type  of 
Surface 

Mesh 

Number  of 
Surface 
Faces 

Number  of 
Prismatic 
Layers 

Radial 

Stretching 

Factor 

Number  of 
Surface 
Nodes 

Total 
Number 
of  Nodes 

“Structured” 

1536 

65 

i.i 

864 

57024 

“Unstructured” 

1604 

65 

i.i 

890 

58740 

Figure  6.7:  Flow  pat  Mines  around  a  cylinder  at  Re  =  40. 


The  pressure  coefficient  on  the  wall  of  the  cylinder  along  the  circumferential 
direction  for  both  grids  is  shown  in  Figure  6.8.  The  pressure  coefficient  is  defined 
in  equation  (6.1). 


P~  Poo 


'v  —  i 


\PooUl 


(6.1) 


The  angle  theta  on  the  abscissa  is  defined  as  the  angle  from  the  forward  stagnation 
point  to  the  aft  stagnation  point  on  the  cylinder.  Numerical  data  are  compared 
with  experimental  data  by  Grove  [46].  The  figure  shows  the  solution  obtained  from 
the  “unstructured”  mesh  compares  well  with  the  experimental  data  and  the  solution 
from  the  “structured”  mesh. 

Nondimension al  vorticity  on  the  wall  of  the  cylinder  is  shown  in  Figure  6.9. 
The  nondimensional  vorticity  is  defined  in  equation  (6.2)  where  D  is  the  diameter 
of  the  cylinder  and  uy  is  the  axial  component  of  the  vorticity. 


.  .  wvD  .  . 

vorticity  =  -j—  (6.2) 

Numerical  results  are  compared  with  the  numerical  data  by  Fornberg  [40].  Form 
berg’s  numerical  data  were  obtained  using  the  stream  function- vorticity  method  with 
130  points  in  the  circumferential  direction  and  114  points  in  the  radial  direction. 
The  figure  shows  the  solution  obtained  from  the  “unstructured”  mesh  compares  rea¬ 
sonably  well  with  the  data  from  the  literature  and  the  solution  from  the  “structured” 
mesh. 

Additional  comparisons  of  the  present  work  on  the  “unstructured”  grid  with 
numerical  and  experimental  results  found  in  the  literature  [20,  36,  97,  131,  133]  are 

shown  in  Table  6.3.  The  separation  angle  is  defined  as  the  angle  from  the  rear 

stagnation  line  at  which  the  attached  boundary  layer  separates  from  the  cylinder. 
The  other  parameters  are  the  nondimensional  length  of  the  bubble,  and  the  drag 
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Figure  6.8:  Pressure  coefficient  distribution  on  cylinder  surface  for  Re=40. 
o  Experiment  by  Grove  [46] 

-  -  Present  work  using  “structured”  grid 

—  Present  work  using  “unstructured”  grid 


Figure  6.9:  Vorticity  distribution  on  cylinder  surface  for  Re=40. 
o  Numerical  results  by  Fornberg  [40] 

-  -  Present  work  using  “structured”  grid 
—  Present  work  using  “unstructured”  grid 
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coefficient  which  is  defined  as  follows: 


Drag 

\PooUlA 


(6.3) 


where  A  is  the  frontal  area  of  the  cylinder  (diameter  times  length).  The  drag 
is  calculated  by  integrating  the  pressure  and  shear  stresses  acting  on  the  faces  of 
the  cells  on  the  cylinder  surface.  The  convergence  of  the  drag  coefficient  and  the 
maximum  residual  of  the  x-momentum  equation  are  shown  in  Figures  6.10  and 
6.11  for  the  case  using  the  “unstructured”  grid.  From  these  results,  the  following 
conclusions  can  be  drawn: 


1.  A  grid  composed  of  prismatic  cells  can  yield  accurate  results  when  compared 
to  experimental  data  and  other  numerical  data  from  the  literature. 

2.  The  pressure  correction  scheme  converges  well  and  yields  stable  solutions  with 
fully  unstructured  prisms. 


3.  The  use  of  a  non-staggered  grid  with  artificial  dissipation  produces  smooth 
pressure  and  vorticity  fields. 

4.  Fully  unstructured  prisms  yield  the  same  accuracy  as  “structured”  prismatic 
meshes. 
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Table  6.3:  Comparisons  of  results  for  circular  cylinder  flow  at  Re=40. 

Separation  Bubble  Length/  Drag 
Angle  (deg)  Cylinder  Diameter  Coefficient 

Present  Work  53  ~  2.0  1.49 

Other  Computations  50  -  53  1.8  -  2.5  1.44  -  1.70 

and  Experiments 


6.3  Impulsive  Start  of  a  Cylinder  using  a  Prismatic 


Mesh 

Impulsive  start  of  a  circular  cylinder  which  attains  a  steady  speed  corresponding  to 
a  Reynolds  number  of  40  is  simulated  in  order  to  examine  the  time  accuracy  of  the 
solver  with  a  prismatic  mesh.  Grids  of  varying  fineness  were  tested  until  the  results 
were  no  longer  mesh  dependent.  The  grid  parameters  used  for  this  case  are  shown 
in  Table  6.4.  The  surface  mesh  is  the  same  fully-unstructured  surface  mesh  used  in 
the  previous  section.  Symmetry  planes  were  defined  at  both  ends  of  the  cylinder. 
The  artificial  dissipation  factor,  a 4,  was  equal  to  10“4  and  the  CFL  factor  was  0.1. 

Initial  conditions  of  the  impulsive  start  simulation  is  potential  flow  [43].  As 
the  simulation  begins,  the  nodes  on  the  aft  stagnation  line  of  the  cylinder  are  mon¬ 
itored  and  the  node  closest  to  the  cylinder  to  register  a  positive  velocity  defines 
the  outermost  position  of  the  separation  bubble.  The  beginning  of  the  separation 
bubble  is  defined  at  the  rear  of  the  cylinder  surface.  The  separation  bubble  length, 
L,  is  nondimensionalized  by  the  cylinder  diameter,  D,  and  the  time  is  the  nondi- 
mensional  time,  tU^/D .  The  growth  of  the  separation  bubble  in  time  is  shown  in 
Figure  6.12.  Experimental  results  from  Coutanceau  [37]  and  Honji  [53]  are  shown 
for  comparison. 

Results  from  the  present  work  match  the  unsteady  portion  of  the  data,  how¬ 
ever,  it  slightly  underpredicts  the  length  of  the  steady  state  value  of  the  separation 
bubble.  Experimental  results  also  show  a  significant  variation  in  the  steady  state 
value  of  the  separation  bubble  size,  and  other  numerical  data  [37]  and  Table  6.3 
show  more  variation.  Other  numerical  results  [33]  report  a  slight  decrease  in  the 
steady  state  value  of  the  separation  bubble  as  the  grid  is  refined.  The  present  work 
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found  a  similar  trend. 


From  these  results,  the  following  conclusion  can  be  drawn:  for  this  impulsive 
start,  the  scheme  is  time  accurate  when  compared  to  experimental  data. 


Table  6.4:  Grid  parameters  used  for  unsteady  circular  cylinder  flow  at  Re=40. 


Number  of 
Surface 
Faces 

Number  of 

Prismatic 

Layers 

Radial 

Stretching 

Factor 

Grid 
Spacing 
at  Wall 

Total 
Number 
of  Nodes 

Radius  of 
Far  Field 
Boundary 

1604 

in 

1.045 

0.01 

99680 

~  20 

95 


Figure  6.12:  Growth  in  time  of  the  separation  bubble  for  an  impulsive  start  of  a 
cylinder  Re  =  40. 

A  Experiment  by  Coutanceau  [37]  o  Experiment  by  Honji  [53] 

—  Present  Work 
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6.4  Steady  3-D  Flow  in  a  Driven  Cubic  Cavity  with  a 
Prismatic  Mesh 

Steady  flow  within  a  driven  cubic  cavity  is  used  to  demonstrate  the  three-dimensional 
capability  of  the  solver  on  a  prismatic  mesh.  A  cubic  cavity  is  modeled  with  sides 
of  length  one.  The  lid  (at  z  =  1.0)  is  moving  with  a  steady  velocity  in  the  positive 
x  direction.  Since  the  flow  is  symmetric  about  the  y=0.5  plane,  only  half  the  cubic 
cavity  is  actually  modeled  with  one  plane  being  a  symmetry  plane.  The  flow  at 
Reynolds  number  of  400,  based  on  the  length  of  the  sides  of  the  cubic  cavity,  is 
modeled.  No  artificial  dissipation  was  used  in  this  case  and  the  CFL  factor  was 
0.1.  Grids  of  varying  fineness  were  tested  until  the  results  were  no  longer  mesh 
dependent.  The  grid  parameters  are  shown  in  Table  6.5. 

Two  velocity  contours  of  the  flow  on  x-z  planes  at  y— 0.015  and  y=0.5  are 
shown  in  Figure  6.13  to  give  a  qualitative  view  of  the  three  dimensionality  of  the 
flow.  The  cavity  lid  is  the  top  surface  and  is  moving  obliquely  from  right  to  left. 
The  grid  mesh  on  the  bottom  wall  is  shown,  with  the  clustered  points  near  walls. 
The  grid  mesh  on  the  x=0  side  wall  and  the  y=0  side  wall  are  also  shown.  Velocity 
contours  in  three  y-z  planes  at  x=0.1,  0.5,  and  0.9  are  shown  in  Figure  6.14.  The 


Table  6.5:  Grid  parameters  used  for  the  cubic  cavity  at  Re  —  400. 


Number  of 
Faces 
on  Wall 

Total 
Number 
of  Nodes 

Number 

of 

Layers 

Smallest 

Grid 

Step  Size 

Largest 

Grid 

Step  Size 

2500 

67626 

50 

0.01 

0.025 
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velocity  contours  do  not  show  the  presence  of  small  vortices  in  the  corners  of  the 
cavity,  however,  velocity  vector  plots,  which  are  not  shown,  indicate  the  vortices  are 
present. 

Two  velocity  profiles  on  the  symmetry  plane  of  the  cubic  cavity  are  plotted 
for  comparison  of  the  solution  accuracy  with  numerical  data  by  Babu  [10].  Figure 
6.15  shows  the  u  velocity  versus  the  z  coordinate  at  x  =  0.5,  and  Figure  6.16  shows 
the  w  velocity  versus  the  x  coordinate  at  z  =  0.5.  Both  plots  show  the  present  work 
compares  well  with  the  numerical  data  from  the  literature.  The  pressure  boundary 
conditions  for  this  cubic  cavity  flow  required  special  consideration.  The  solid  walls 
of  the  cavity  required  boundary  conditions  of 

d<f)  At  d2un 

8n  Re  On 2  V  ' 

where  n  represents  the  normal  direction  to  the  specific  wall.  Derivatives  were  either 
forward  or  backward  differenced,  depending  on  the  location  of  the  wall.  From  these 
results,  the  following  conclusions  can  be  drawn: 

1.  The  simulated  flow  within  a  three-dimensional  geometry,  using  a  grid  com¬ 
posed  of  prismatic  cells,  contains  the  physical  characteristics  of  the  expected 
flow. 

2.  The  solver  yields  accurate  results  when  compared  to  other  numerical  data  from 
the  literature. 
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Figure  6.15:  U  velocity  profile  versus  z  on  the  symmetry  plane  at  x  =  0.5  for  the 
cubic  cavity  at  Re  =  400. 

o  Numerical  result  by  Babu  [10]  - Present  work 
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Figure  6.16:  W  velocity  profile  versus  x  on  the  symmetry  plane  at  z  =  0.5  for  the 
cubic  cavity  at  Re  =  400. 

o  Numerical  result  by  Babu  [10]  - Present  work 
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6.5  Steady  3-D  Flow  Around  a  Sphere  using  a  Prismatic 


Mesh 

Steady  flow  around  a  sphere  is  used  to  demonstrate  the  three-dimensional  capability 
of  the  solver  on  a  prismatic  mesh.  The  spherical  flow  is  calculated  for  several 
Reynolds  numbers,  based  on  the  diameter  of  the  sphere.  At  very  low  Reynolds 
numbers,  the  viscous  terms  dominate  and  the  convective  terms  are  negligible  (Stokes 
flow).  As  the  Reynolds  number  increases,  convective  terms  become  significant  and 
Oseen  developed  a  correction  to  the  equations  for  Stokes  flow.  When  the  Reynolds 
number  exceeds  1.0,  both  Stokes  and  Oseen’s  flow  equations  are  invalid,  and  the 
resulting  flow  is  steady  and  attached.  An  attached  axisymmetric  ring  eddy  appears 
when  the  Reynolds  number  is  approximately  20  [139].  The  ring  eddy  remains  steady 
and  attached  until  the  Reynolds  number  reaches  130.  A  slight  oscillation  in  the  wake 
[92,  117]  develops  at  this  Reynolds  number.  Experiments  [117]  have  shown  periodic 
hairpin  vortex  shedding  begins  at  approximately  a  Reynolds  number  of  300.  These 
vortices  produce  an  unsteady  laminar  wake.  When  the  Reynolds  number  is  increased 
to  800,  the  wake  becomes  turbulent  and  the  structure  of  the  vortices  become  unclear. 

At  low  Reynolds  numbers,  below  300,  the  flow  over  a  sphere  is  basically 
steady  and  symmetric.  Therefore,  the  flow  simulations  at  these  Reynolds  numbers 
used  a  half  sphere  on  a  symmetry  plane.  Grids  of  varying  fineness  were  tested  until 
the  results  were  no  longer  mesh  dependent.  Grid  parameters  used  for  the  sphere 
calculations  are  shown  in  Table  6.6.  All  the  grids  had  1346  triangles  on  the  surface 
mesh  and  used  a  radial  stretching  factor  of  1.1.  The  convergence  criterion  for  the 
maximum  residual  in  the  x-momentum  equation  was  10-5.  The  artificial  dissipation 
factor,  <t4,  was  equal  to  10~4  and  the  CFL  factor  was  0.1.  The  prismatic  mesh  for 
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Table  6.6:  Prismatic  grid  parameters  used  for  flow  around  a  sphere. 


Reynolds 

Number 

Number  of 
Prismatic 
Layers 

Grid 
Spacing 
at  Wall 

Total 
Number 
of  Nodes 

Radius  of 
Far  Field 
Boundary 

0.1 

57 

0.02 

40890 

~  20 

1.0 

65 

0.01 

46530 

~  21 

10.0 

65 

46530 

~  21 

40.0 

65 

46530 

~  21 

100.0 

65 

46530 

~  21 

200.0 

65 

0.01 

46530 

~  21 

the  Reynolds  number  equal  to  10.0  case  is  shown  in  Figure  6.17. 

Classic  experimental  drag  coefficient  versus  Reynolds  number  data  from 
Schlichting  [120],  as  well  as  other  experimental  data  with  an  error  range  [115],  are 
used  for  comparison,  as  well  as  the  analytic  results  from  Stokes  flow  around  a  sphere 
and  Oseen’s  equation  for  a  sphere  [97,  120].  The  drag  coefficient  in  the  present  work 
is  defined  by  the  following  equation: 


Drag 

%PooUZ,A 


(6.5) 


where  A  is  the  frontal  area  of  the  sphere  (it  r 2)  where  r  is  the  radius  of  the  sphere. 


Figure  6.17:  Prismatic  surface  and  symmetry  plane  mesh  for  sphere  at  Re 


The  drag  coefficient  for  Stokes  flow  is  given  as  follows,  while  the  coefficient  corre¬ 
sponding  to  Oseen’s  solution  is  given  by  the  following  equation. 


CdStokes  =  (6.6) 

_  24.0  /  3  „  \  ,  , 

C dO seen  -  Re  ^  +  jg^J  (6.7) 

Results  of  the  drag  coefficient  comparisons  are  shown  in  Figure  6.18.  Numerical 
results  agree  quite  well  with  the  theoretical  results  and  fall  within  the  scatter  of  the 
experimental  results. 

The  wall  pressure  coefficient  and  vorticity  are  plotted  in  Figures  6.19  and 
6.20,  respectively,  for  the  Re=0.01  case.  The  pressure  coefficient  shows  an  antisym¬ 
metric  profile  which  is  characteristic  of  Stokes  flow.  The  vorticity  shows  a  symmetric 
profile  which  is  also  characteristic  of  Stokes  flow,  ie.,  no  separation  has  occurred. 
The  wall  pressure  coefficient  and  vorticity  are  plotted  in  Figures  6.21  and  6.22,  re¬ 
spectively,  for  the  Re— 10,  40,  100,  and  200  cases.  Both  figures  show  separation 
occurring  on  the  aft  side  of  the  sphere  when  the  Reynolds  number  is  greater  than 
40.  The  pathlines  around  the  sphere  are  illustrated  for  the  Reynolds  number  equal 
to  200  case  in  Figures  6.23  and  6.24.  Two  orthogonal  plane  cuts  are  shown  to  illus¬ 
trate  the  presence  of  the  ring  eddy  on  the  rear  of  the  sphere.  From  these  results, 
the  following  conclusions  can  be  drawn: 

1.  The  incompressible  solver  yields  accurate  solutions,  for  flow  over  a  three- 
dimensional  body,  when  compared  to  experimental  data. 

2.  The  solver  yields  expected  results  over  a  range  of  Reynolds  numbers  which 
includes  Stokes  flow  and  separated  flows. 
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Figure  6.18:  Drag  coefficient  of  a  sphere  versus  Reynolds  number. 
Tic  at  lower  right  represents  error  range  of  experimental  data. 

- Stokes  flow 

- Oseen’s  equation 

- Experimental  data  from  Schlichting  [120]  and  Roos  [115] 

o  Present  work  with  hemisphere  grid 
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Figure  6.19:  Pressure  coefficient  distribution  on  the  surface  of  the  sphere  for  Re=0.1 
(Stokes  flow). 


Figure  6.20:  Vorticity  distribution  on  the  surface  of  the  sphere  for  Re=0.1  (Stokes 
flow). 
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e  surface  of  the  sphere  for  low  Re. 


Figure  6.24:  Pathlines  on  the  x-y  plane  near  a  sphere  at  Re  =  200. 
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6.6  Evaluation  of  Prismatic  Mesh  Adaptation  Applied 


to  Cubic  Cavity  Flow 

Prismatic  adaptation  will  be  evaluated  by  examining  the  accuracy  of  the  solution 
with  local  embedding  compared  to  the  solution  from  an  unembedded  coarse  grid. 
The  solution  obtained  with  local  embedding  should  be  the  same  as  the  solution 
obtained  by  the  corresponding  globally  embedded  mesh.  Robustness  of  the  scheme 
will  also  be  examined. 

6.6.1  Accuracy  of  Prismatic  Adaptation  on  Cubic  Cavity  Flow 

Steady  flow  in  a  cubic  cavity  at  a  Reynolds  number  of  400  is  calculated  using  coarse, 
adapted,  and  globally  adapted  prismatic  grids.  The  cubic  cavity  is  modeled  with 
sides  of  length  one.  The  lid  (at  z  =  1.0)  is  moving  with  a  steady  velocity  in  the 
positive  x  direction.  Since  the  flow  is  symmetric  about  the  y=0.5  plane,  only  half 
the  cubic  cavity  is  actually  modeled  with  one  plane  being  a  symmetry  plane.  All 
three  grids  have  the  same  point  distribution  in  the  z  direction  as  shown  in  Table 
6.5.  Only  the  surface  triangulation  is  different  in  the  three  cases.  The  adapted 
grid  was  obtained  from  the  coarse  grid  using  the  solution  after  5000  iterations.  The 
CFL  factor  for  all  three  cases  was  0.1,  and  no  artificial  dissipation  was  used.  The 
grid  parameters  for  the  three  grids  are  shown  in  Table  6.7.  The  coarse,  adapted, 
and  globally  adapted  grids  are  shown  in  Figures  6.25,  6.26,  and  6.27,  respectively. 
In  Figure  6.26,  directional  embedding  of  the  prisms  is  focused  near  the  side-walls 
of  the  cavity.  This  is  due  to  the  flow  decelerating  near  the  wall  boundaries.  The 
boundary  which  is  not  adapted  is  a  symmetry  plane.  The  solid  walls  of  the  cavity 
had  boundary  conditions  of  =  0.0  where  n  represents  the  normal  direction  to  the 
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Figure  6.25:  Coarse  grid  surface  triangulation  for  cubic  cavity  at  Re  =  400. 


specific  wall. 

Two  velocity  profiles  on  the  symmetry  plane  of  the  cubic  cavity  are  plotted 
for  comparison  of  the  solution  accuracy.  Figure  6.28  shows  the  u  velocity  versus  the 
z  coordinate  at  x  =  0.5,  and  Figure  6.29  shows  the  w  velocity  versus  the  x  coordinate 
at  z  =  0.5.  There  is  no  significant  difference  between  the  three  solutions  in  Figure 
6.28,  but  the  coarse  grid  gives  a  very  poor  solution  for  the  w  velocity  in  Figure  6.29. 


6.6.2  Robustness  and  CPU  Time  Savings  of  Prismatic  Adaptation 

Robustness  of  the  prismatic  adapter  is  demonstrated  in  Figure  6.30.  The  maximum 
residual  in  the  x-momentum  equation  is  plotted  versus  time  for  the  adapted  and 
globally  adapted  cases.  The  large  spike  in  Figure  6.30  occurs  due  to  the  linearly 
interpolated  solution  from  the  initial  coarse  grid  nodes  to  the  newly  created  ones 
due  to  adaptation.  The  flow  solver  converges  relatively  quickly  after  the  restart  to 
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Figure  6.25:  Coarse  grid  surface  triangulation  for  cubic  cavity  at  Re  =  400. 


specific  wall. 

Two  velocity  profiles  on  the  symmetry  plane  of  the  cubic  cavity  are  plotted 
for  comparison  of  the  solution  accuracy.  Figure  6.28  shows  the  u  velocity  versus  the 
z  coordinate  at  x  =  0.5,  and  Figure  6.29  shows  the  w  velocity  versus  the  x  coordinate 
at  z  =  0.5.  There  is  no  significant  difference  between  the  three  solutions  in  Figure 
6.28,  but  the  coarse  grid  gives  a  very  poor  solution  for  the  w  velocity  in  Figure  6.29. 


6.6.2  Robustness  and  CPU  Time  Savings  of  Prismatic  Adaptation 

Robustness  of  the  prismatic  adapter  is  demonstrated  in  Figure  6.30.  The  maximum 
residual  in  the  x-momentum  equation  is  plotted  versus  time  for  the  adapted  and 
globally  adapted  cases.  The  large  spike  in  Figure  6.30  occurs  due  to  the  linearly 
interpolated  solution  from  the  initial  coarse  grid  nodes  to  the  newly  created  ones 
due  to  adaptation.  The  flow  solver  converges  relatively  quickly  after  the  restart  to 
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Figure  6.28:  Adaptation  comparison  of  the  u  velocity  profile  versus  z  on  the  sym¬ 
metry  plane  at  x  =  0.5  for  the  cubic  cavity  at  Re  =  400. 

- Coarse  grid  - Locally  adapted  grid 

- Globally  adapted  grid 


Figure  6.29:  Adaptation  comparison  of  the  w  velocity  profile  versus  x  on  the  sym¬ 
metry  plane  at  z  =  0.5  for  the  cubic  cavity  at  Re  =  400. 

- Coarse  grid  - Locally  adapted  grid 

- Globally  adapted  grid 
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Table  6.7:  Surface  triangulation  used  for  adaptation  accuracy  comparison  for  the 
cubic  cavity  at  Re  =  400. 


Type  of  Grid 

Number  of  Surface  Faces 

Total  Number  of  Nodes 

Coarse 

994 

27693 

Locally  Adapted 

3068 

82518 

Globally  Adapted 

3976 

106029 

residual  levels  on  the  order  of  10“5. 

Another  criterion  used  to  evaluate  an  algorithm  is  efficiency.  Computing 
time  to  reach  steady  state  using  the  locally  adapted  grid  as  well  as  the  globally 
adapted  grid  are  compared.  The  CPU  time  for  the  adapted  grid  consists  of  the 
CPU  time  for  the  initial  5000  iterations  with  the  coarse  grid,  and  the  CPU  time  for 
the  adapted  grid  for  the  final  5000  iterations.  The  globally  adapted  fine  grid  CPU 
time  is  for  10000  iterations  which  were  required  to  reach  steady  state  conditions. 
The  maximum  residuals  in  the  x-momentum  equation  were  on  the  order  of  10~6. 
The  CPU  time  for  the  locally  adapted  grid  and  the  globally  adapted  grid  are  shown 
in  Table  6.8.  The  CPU  time  for  the  locally  adapted  grid  case  will  vary  depending 
on  the  time  the  coarse  grid  is  adapted  and  the  time  the  simulation  is  terminated.  In 
this  test  case,  there  is  more  than  50%  CPU  time  savings  using  the  locally  adapted 
grid  for  the  same  solution  accuracy, 

From  these  results,  the  following  conclusions  can  be  drawn: 
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Figure  6.30:  Demonstration  of  robustness  of  prismatic  adapter,  convergence  of  the 
maximum  residual  of  the  x-momentum  equation  for  the  cubic  cavity  at  Re  =  400. 

-  -  -  Locally  adapted  grid  - Globally  adapted  grid 
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Table  6.8:  Adaptation  CPU  time  comparison  for  the  cubic  cavity  at  Re  =  400. 


Type  of  Grid 

Number  of  Cells 

CPU  Time  (sec) 

Locally  Adapted 

153400 

31300 

Globally  Fine 

198800 

64100 

1.  The  solver  with  locally  adapted  prismatic  cells  yields  a  stable  solution. 

2.  The  locally  adapted  prismatic  mesh  yields  the  same  accuracy  with  reduced 
computing  resources  compared  to  a  globally  adapted  prismatic  mesh. 
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6.7  Equation  Adaptation  with  a  Hybrid  Mesh  for  Flow 


Around  a  Sphere 

Equation  adaptation  of  the  solver  is  evaluated  for  the  case  of  flow  over  a  sphere 
at  Re  =  40  using  a  hybrid  mesh.  Two  cases  are  considered  here.  The  first  case 
solved  the  Euler  equations  in  the  tetrahedral  region  (“inviscid”  tetrahedra),  while 
the  second  case  solved  the  Navier-Stokes  equations  in  the  tetrahedral  region  (“ 
viscous”  tetrahedra).  The  full  Navier-Stokes  equations  are  employed  over  the  prism 
region  for  both  cases.  The  grid  parameters  used  for  this  study  are  shown  in  Table 
6.9.  Both  simulations  terminated  after  6000  iterations  with  the  maximum  residuals 
in  the  x-momentum  equation  being  on  the  order  of  10”4.  The  CFL  factor  is  0.1,  and 
the  momentum  artificial  dissipation  parameter,  <74,  is  10“4.  The  pressure  coefficient 
and  vorticity  are  examined  along  the  symmetry  plane  (equator)  of  the  sphere  in 
Figures  6.31  and  6.32.  The  angle  theta  is  defined  as  zero  at  the  front  stagnation 
point  and  180  degrees  at  the  rear  stagnation  point  of  the  sphere.  It  is  observed  that 
the  solutions  corresponding  to  the  two  cases  are  nearly  identical.  The  CPU  times 
for  the  two  simulations  are  shown  in  Table  6.10.  The  initial  case  with  “inviscid 
tetrahedra”  requires  about  60  %  of  the  time  required  by  the  case  that  does  not  use 
equation  adaptation.  It  should  be  noted  that  the  hybrid  (prismatic/tetrahedral) 
mesh  makes  equation  adaptation  simple  to  implement. 

The  effect  of  the  interface  between  the  prisms  and  the  tetrahedra  on  the 
solution  is  examined  next.  Flow  around  the  sphere  at  Re  =  100  is  considered.  The 
prismatic  region  was  reduced  so  that  the  separation  bubble  would  no  longer  lie  solely 
in  the  “viscous”  prism  region.  Figure  6.33  shows  the  pathlines  around  the  sphere. 
It  is  observed  that  the  pathlines  are  smooth  across  the  interface.  Figures  6.34  and 
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Table  6.9:  Grid  parameters  \ised  for  evaluation  of  equation  adaptation  for  sphere 
flow  at  Re=40. 


Number 

of 

Surface 

Faces 

Number 

of 

Prismatic 

Layers 

Radial 

Stretching 

Factor 

Grid 

Spacing 

at 

Wall 

Number 

of 

Prism 

Nodes 

Total 

Number 

of 

Nodes 

Number 

of 

Tetrahedral 

Cells 

1346 

27 

i.i 

0.01 

19740 

49195 

160054 

6.35  illustrate  the  velocity  and  pressure  coefficient  contours  on  the  symmetry  plane. 

The  contours  are  relatively  smooth  across  the  interface. 

From  these  results,  the  following  conclusions  can  be  drawn: 

1.  The  solver  with  equation  adaptation  and  a  hybrid  (prism/tetrahedra)  mesh 
yields  stable  solutions. 

2.  The  solver  with  equation  adaptation  yields  nearly  the  same  results  with  re¬ 
duced  computing  resources  compared  to  the  case  with  solves  the  Navier- Stokes 
equations  over  the  entire  mesh. 

3.  The  change  of  topology  of  the  grid  (prism/tetrahedra  interface)  has  little  effect 
on  the  solution. 
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Table  6.10:  Equation  adaptation  CPU  time  comparison  for  the  sphere  at  Re  =  40. 


Type  of  Tetrahedra 

CPU  Time  (sec) 

“Inviscid” 

29800 

“Viscous” 

49700 

Figure  6.31:  Evaluation  of  equation  adaptation  with  hybrid  mesh.  Pressure  coeffi¬ 
cient  distributions  on  the  surface  of  the  sphere  for  Re=40. 

- “Inviscid”  Tetrahedra 

- “Viscous”  Tetrahedra 
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Figure  6.32:  Evaluation  of  equation  adaptation  with  hybrid  mesh.  Vorticity  distri¬ 
butions  on  the  surface  of  the  sphere  for  Re=40. 

-  “Inviscid”  Tetrahedra 

- “Viscous”  Tetrahedra 
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Figure  6.34:  Effect  of  prism/tetrahedra  interface  on  the  solution.  Velocity  contours 
near  a  sphere  for  Re=100  are  relatively  unaffected  by  the  change  in  topology  of  the 
grid. 
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Figure  6.35:  Effect  of  prism/tetrahedra  interface  on  the  solution.  Pressure  contours 
near  a  sphere  for  Re=100  are  relatively  unaffected  by  the  change  in  topology  of  the 


6.8  Hybrid  Mesh  Adaptation  Applied  to  Flow  Around 


a  Sphere 

Steady  flow  over  a  sphere  is  used  to  demonstrate  the  advantages  of  hybrid  mesh 
adaptation.  The  hybrid  mesh  consists  of  “viscous”  prisms  near  the  surface  of  the 
sphere,  and  “inviscid”  tetrahedra  filling  the  rest  of  the  computational  domain.  The 
robustness  of  the  scheme  will  also  be  presented. 

6.8.1  Accuracy  of  Hybrid  Adaptation  on  Sphere  Flow 

Steady  flow  over  a  sphere  at  a  Reynolds  number  of  100  is  used  to  demonstrate  the 
accuracy  of  hybrid  mesh  adaptation.  The  parameters  of  the  coarse,  locally  adapted, 
and  globally  fine  grids  used  for  this  study  are  shown  in  Table  6.11.  The  globally 
fine  mesh  was  created  by  globally  refining  the  prismatic  region  of  the  coarse  grid. 
Refinement  occurred  on  portions  of  the  tetrahedral  mesh  due  to  the  refinement  of 
the  prisms.  In  the  prism  region,  all  three  grids  had  36  layers,  with  the  grid  step 
size  at  the  wall  being  0.01,  and  a  stretching  factor  of  1.1.  The  artificial  dissipation 
factor,  04,  was  equal  to  10-4  and  the  CFL  factor  was  0.1. 

Wall  pressure  coefficient  and  vorticity  are  shown  in  Figures  6.36  and  6.37. 
The  pressure  coefficient  solution  for  the  adapted  mesh  and  the  fine  mesh  are  nearly 
identical.  However,  the  vorticity  distribution  corresponding  to  the  adapted  and  fine 
meshes  have  some  differences.  The  two  solutions  agree  over  the  rear  portion  of  the 
sphere,  while  differences  are  observed  near  6  =  60  degrees.  An  examination  of  the 
adapted  mesh  shows  the  rear  of  the  sphere  is  adapted,  while  the  region  near  ^  =  60 
degrees  was  not  adapted.  At  low  Reynolds  numbers,  the  sphere  actually  requires  a 
globally  fine  mesh.  However,  this  comparison  of  adapted  and  fine  grids  still  shows 
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the  benefit  of  adaptation. 

The  adapted  grid  is  shown  in  Figure  6.39  along  with  pathlines  around  the 
sphere.  The  figure  shows  the  tessellation  on  the  wall  surface  and  on  the  symmetry 
plane.  It  is  clearly  seen  that  embedding  in  the  tetrahedral  region  is  focused  at  the 
rear  of  the  sphere.  The  prismatic  region  is  also  directionally  refined  near  the  up¬ 
stream  and  downstream  sections  of  the  body.  This  is  due  to  the  flow  accelerating 
from  the  upstream  stagnation  point  and  the  flow  downstream  separating  causing 
significant  flow  gradients  in  the  lateral  directions  which  are  detected  by  the  direc¬ 
tional  adaptive  algorithm.  Pathlines  show  the  ring  eddy  formed  at  the  rear  of  the 
sphere,  and  the  cells  in  the  tetrahedral  region  of  the  hybrid  mesh.  Figure  6.40  shows 
a  different  view  of  the  adapted  hybrid  mesh  along  with  velocity  contours.  Three 
surfaces  are  shown.  The  first  is  the  sphere  surface,  the  second  is  the  symmetry 
plane,  and  the  third  is  a  plane  cutting  through  the  interior  of  the  grid,  normal  to 
the  symmetry  plane.  Despite  the  drastic  changes  in  topology  at  the  interface,  the 
contours  lines  across  it  are  relatively  smooth.  Figure  6.41  shows  the  pressure  coef¬ 
ficient  contours  on  the  symmetry  plane.  The  pressure  field  is  relatively  smooth  as 
well  as  the  contour  lines  across  the  prismatic/tetrahedral  interface. 

6.8.2  Robustness  and  CPU  Time  Savings  of  Hybrid  Adaptation 

Robustness  of  the  hybrid  adapter  is  demonstrated  in  Figure  6.38.  The  maximum 
residual  of  the  x-momentum  equation  is  plotted  versus  time  for  the  locally  adapted 
and  globally  adapted  cases.  The  large  spike  in  Figure  6.38  is  due  to  the  linearly 
interpolated  solution  from  the  initial  coarse  grid  nodes  to  the  newly  created  ones 
due  to  adaptation.  The  flow  solver  converges  relatively  quickly  after  the  restart  to 
residual  levels  on  the  order  of  10~4. 


125 


Another  criteria  used  to  evaluate  an  algorithm  is  efficiency.  The  CPU  time 
to  reach  steady  state  for  the  locally  adapted  grid  and  the  globally  adapted  grid  are 
compared  in  Table  6.12.  The  CPU  time  for  the  locally  adapted  grid  case  will  vary 
depending  on  the  time  the  coarse  grid  is  adapted  and  the  time  the  simulation  is 
terminated.  For  this  case,  there  is  more  than  40  %  CPU  time  savings  using  the 
locally  adapted  grid. 

From  these  results,  the  following  conclusions  can  be  drawn: 

1.  The  solver  with  hybrid  adaptation  yields  stable  solutions. 

2.  The  locally  adapted  hybrid  mesh  yields  approximately  the  same  results  with 
reduced  computing  resources  compared  to  a  fine  mesh. 

3.  The  incompressible  pressure  correction  method  with  non-staggered  grids  yields 
smooth  pressure  and  velocity  fields. 

4.  The  velocity  field  and  flow  pathlines  are  unaffected  by  the  change  in  grid 
topology  at  the  prism  and  tetrahedra  interface. 
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Table  6.11:  Grid  parameters  used  for  hybrid  mesh  adaptation  for  sphere  flow  at 
Re=100. 


Type 

Number  of 

Total 

Number 

Number  of 

of 

Surface 

Number 

of  Prism 

Tetrahedral 

Grid 

Faces 

of  Nodes 

Cells 

Cells 

Coarse 

332 

13372 

11952 

35063 

Locally  Adapted 

1095 

31981 

39420 

59401 

Fine 

1328 

36495 

47808 

59796 

Theta  (degrees) 


Figure  6.36:  Pressure  coefficient  distribution  on  the  surface  of  a  sphere  with  a  hybrid 
adapted  mesh  for  Re=100. 

-  -  -  Coarse  grid  - Locally  adapted  grid 

- Fine  grid  (globally  adapted  prism  region) 
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Figure  6.37:  Vorticity  distribution  on  the  surface  of  a  sphere  with  a  hybrid  adapted 
mesh  for  Re=100. 

- Coarse  grid  - Locally  adapted  grid 

- Fine  grid  (globally  adapted  prism  region) 


Time 


Figure  6.38:  Demonstration  of  robustness  of  hybrid  adapter,  maximum  residual  in 
the  x-momentum  equation  versus  time  for  the  sphere  at  Re  =  100. 

- Locally  adapted  grid  - Globally  adapted  grid 
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Figure  6.39:  Adapted  hybrid  mesh  and  flow  pathlines  for  sphere  at  Re=100.  A 
view  of  the  tessellation  on  the  wall  surface  and  symmetry  plane.  The  hybrid  grid  is 
embedded  isotropically  in  the  tetrahedral  region  and  directionally  in  the  prismatic 
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Figure  6.40:  Adapted  hybrid  mesh  and  corresponding  flow  velocity  contours  for 
sphere  at  Re=100.  A  view  of  the  tessellation  on  the  wall  surface,  symmetry  plane 
and  an  interior  ecpiatorial  plane.  The  hybrid  grid  is  embedded  isotropically  in  the 
tetrahedral  region  and  directionally  in  the  prismatic  region. 
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Figure  6.41:  Adapted  hybrid  mesh  and  pressure  coefficient  contours  for  sphere  at 
Re=100.  A  view  of  the  tessellation  on  the  wall  surface  and  symmetry  plane.  The 
hybrid  grid  is  embedded  isotropically  in  the  tetrahedral  region  and  directionally  in 
the  prismatic  region. 
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Table  6.12:  Adaptation  CPU  comparison  for  the  hybrid  sphere  mesh  at  Re  =  100. 


Type  of  Grid 

Number  of  Nodes 

CPU  Time  (sec) 

Locally  Adapted 

31981 

11510 

Globally  Fine 

36495 

19720 
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6.9  Hybrid  Mesh  Adaptation  for  Flow  Around  Tandem 


Spheres 

Hybrid  mesh  adaptation  is  used  to  simulate  the  steady  flow  around  two  spheres  in 
a  tandem  (in-line)  arrangement.  A  symmetry  plane  is  used  to  model  half  of  the 
flow  domain  at  Re  =  100.  The  grid  parameters  are  shown  in  Table  6.13.  The  grid 
had  a  step  size  at  the  wall  of  0.01  with  a  stretching  factor  of  1.1.  The  surface 
triangulation  and  symmetry  plane  of  the  hybrid  grid  are  shown  in  Figure  6.42.  The 
prismatic  region  shows  directional  adaptation  on  the  fore  and  aft  portions  of  both 
spheres  due  to  the  flow  accelerating  from  the  upstream  stagnation  point  and  the 
flow  downstream  separating  or  decelerating.  Embedding  in  the  tetrahedral  region 
is  focused  at  the  rear  of  the  second  sphere  as  well  as  the  region  between  the  two 
spheres.  The  artificial  dissipation  factor,  <74,  was  equal  to  10~5,  the  CFL  factor  was 
0.05,  and  the  maximum  residual  in  th  x-momentum  equation  was  on  the  order  of 
10"5. 

The  drag  coefficient  on  the  first  sphere  was  1.26,  while  the  drag  coefficient  on 
the  second  sphere  was  0.8.  The  drag  coefficient  for  a  isolated  sphere  is  approximately 
1.1.  The  first  sphere  apparently  “shields”  the  second  sphere  from  the  freestream  flow, 
resulting  in  lower  drag  on  the  second  sphere.  This  shielding  phenomenon  was  seen 
experimentally  [142]  for  tandem  cylinders  in  oscillating  flows.  The  presence  of  the 
second  sphere  inhibits  the  formation  of  the  vortex  aft  of  the  first  sphere,  although  the 
flow  in  this  region  does  appear  to  be  separated.  This  can  be  seen  when  examining 
the  pathlines  shown  in  Figure  6.42.  The  shielding  of  the  second  sphere  by  the  first 
sphere  results  in  the  flow  around  the  second  sphere  to  be  representative  of  lower 
Reynolds  number  flow.  In  fact,  the  flow  past  the  second  sphere  does  not  separate 
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Table  6.13:  Adapted  hybrid  grid  parameters  used  for  tandem  spheres  flow  at 
Re=100. 


Number  of 
Surface 
Faces 

Number 
of  Prism 
Nodes 

3470 

10 

18070 

74309 

34700 

413938 

when  it  reaches  the  aft  portion  of  the  sphere.  A  closeup  view  of  the  region  between 
the  two  spheres  is  shown  in  Figure  6.43.  The  surface  triangulation  of  the  fore  and  aft 
spheres  are  shown  on  the  left  and  right  edges  of  the  figure.  The  quadrilateral  faces 
of  the  prismatic  region  are  shown  as  well  as  the  triangular  faces  of  the  tetrahedral 
cells  on  the  symmetry  plane.  The  velocity  vectors  in  this  region  are  shown  and  the 
flow  is  separated  aft  of  the  first  sphere.  To  simplify  the  figure,  only  velocity  vectors 
with  magnitudes  less  than  0.015  are  shown.  A  different  view  of  the  adapted  hybrid 
mesh,  as  well  as  of  the  flow  velocity  contours  are  shown  in  Figure  6.44. 

From  these  results,  the  following  conclusions  can  be  drawn: 

1.  The  hybrid  mesh  is  suitable  to  cover  multiple  three-dimensional  bodies. 

2.  The  adapted  hybrid  mesh  yields  a  stable  solution  for  this  two-body  problem. 

3.  The  presence  of  the  second  sphere  inhibits  the  formation  of  the  vortex  aft  of 
the  first  sphere. 

4.  The  presence  of  the  first  sphere  “shields”  the  second  sphere  from  the  main 
flow  resulting  in  reduced  drag. 
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Figure  6.42:  Adapted  hybrid  grid  and  flow  pathlines  for  tandem  spheres  at  Re=100. 
A  view  of  the  tessellation  on  the  wall  surface  and  symmetry  plane.  The  hybrid  grid 
is  embedded  isotropically  in  the  tetrahedral  region  and  directionally  in  the  prismatic 


Figure  6.43:  Velocity  vectors  on  the  symmetry  plane  between  tandem  spheres  at 
Re=100.  The  flow  is  separated  aft  of  the  first  sphere.  The  closeup  view  shows  a 
portion  of  the  surface  mesh  of  the  fore  and  aft  spheres,  the  quadrilateral  faces  of  the 
prismatic  cells  on  the  symmetry  plane,  and  the  triangular  faces  of  the  tetrahedral 
cells  on  the  symmetry  plane.  Only  the  velocity  vectors  with  magnitudes  less  than 
0.015  are  shown. 


Figure  6.44:  Adapted  hybrid  grid  and  flow  velocity  contours  for  tandem  spheres  at 
Re=100.  A  view  of  the  tessellation  on  the  wall  surface,  symmetry  plane  and  an  in¬ 
terior  equatorial  plane.  The  hybrid  grid  is  embedded  isotropically  in  the  tetrahedral 
region  and  directionally  in  the  prismatic  region. 
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Chapter  7 


Concluding  Remarks 

7.1  Summary 

The  main  goal  of  the  present  work  was  to  develop  an  adaptive  numerical  algorithm 
for  the  simulations  of  three-dimensional  viscous  incompressible  flows.  The  method 
should  be  stable  at  high  Reynolds  numbers,  accurate,  and  flexible  to  use  with  com¬ 
bined  prismatic  and  tetrahedral  grids. 

A  hybrid  adaptive  solver  for  three-dimensional  incompressible  viscous  flow 
simulations  was  developed  and  validated.  The  hybrid  grid,  comprising  of  prisms  and 
tetrahedra,  proved  to  be  appropriate  for  these  computations.  Equation  adaptation 
resulted  in  appreciable  reduction  in  computing  time.  The  interface  between  the 
prisms  and  the  tetrahedra  did  not  affect  accuracy.  The  semi-structured  nature  of 
the  prismatic  grids  reduced  the  memory  requirements  considerably  as  opposed  to 
using  a  fully  unstructured  grid.  A  dual  grid  adaptation  algorithm  was  implemented 
in  conjunction  with  the  solver.  The  algorithm  enabled  the  solver  to  obtain  accurate 
results  adaptively  using  relatively  fewer  grid  nodes. 
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7.2  Contributions 


The  contributions  of  the  present  work  can  be  categorized  into  three  main  areas. 

These  are  integration  scheme,  unstructured  mesh,  and  adaptation. 

1.  A  time  accurate  finite  volume/finite  element  integration  scheme  for  solving  the 
incompressible  Navier-Stokes  equations  on  hybrid  grids  was  developed.  An 
important  feature  of  the  integration  scheme  is  the  strategy  of  performing  all 
the  operations  on  an  edge- wise  basis  which  makes  it  computationally  efficient, 
especially  for  hybrid  grids  which  consist  of  different  types  of  elements.  This  is 
due  to  the  fact  that  performing  the  operations  on  an  edge-wise  basis  enables 
distinction  between  different  types  of  cells  to  be  ignored  and  hence  makes  the 
scheme  very  general. 

2.  The  solver  has  been  developed  using  non-staggered  grids  rather  than  the  clas¬ 
sical  staggered  grid  approach.  Artificial  dissipation  used  in  conjunction  with 
non-staggered  grids  produced  non-oscillatory  solutions. 

3.  This  is  the  first  incompressible  flow  method  with  prismatic  and  hybrid  grids. 

4.  An  equation  adaptation  scheme  was  coupled  with  the  hybrid  mesh  such  that 
the  Navier-Stokes  equations  are  solved  over  the  prism  region  and  the  simpler 
Euler  equations  are  employed  over  the  tetrahedral  region. 

5.  This  is  the  first  incompressible  flow  method  with  hybrid  adaptation. 

7.3  Conclusions 

The  following  conclusions  can  be  drawn  from  the  present  work. 
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1.  The  incompressible  three-dimensional  Navier-Stokes  solver  is  relatively  accu¬ 
rate  on  stretched  and  skewed  unstructured  grids.  The  scheme  is  also  time 
accurate.  Comparisons  were  made  with  existing  experimental  observations 
where  possible,  and  other  numerical  observations  in  other  cases.  Surface  pres¬ 
sure  and  vorticity  distributions  compared  well  with  the  literature,  as  well  as 
steady  state  sphere  drag  results.  Simulation  of  impulsive  start  of  flow  around 
a  cylinder  demonstrated  time  accuracy. 

2.  The  scheme  is  stable  at  high  Reynolds  numbers.  Stability  at  high  Reynolds 
numbers  is  a  concern  for  all  Navier-Stokes  solvers.  Flow  over  a  flat  plate  at  a 
high  Reynolds  number  was  completed  to  demonstrate  stability  of  the  solver. 

3.  The  upwind-like  artificial  dissipation  model  was  effective  in  eliminating  the 
odd-even  modes  in  the  solution.  Pressure  correction  methods  using  non-staggered 
grids  are  noted  for  odd-even  modes  arising  in  the  solution.  The  fourth  order 
artificial  dissipation  successfully  suppressed  these  modes. 

4.  Use  of  equation  adaptation ,  viscous  prisms  and  inviscid  tetrahedra,  can  yield 
accuracy  that  is  comparable  to  solving  the  Navier-Stokes  equations  over  the 
entire  domain.  Equation  adaptation  was  demonstrated  on  a  hybrid  sphere 
mesh.  Solutions  from  both  grids  were  nearly  identical. 

5.  Use  of  equation  adaptation,  viscous  prisms  and  inviscid  tetrahedra ,  provides 
efficiency  by  yielding  substantial  CPU  time  and  memory  savings  compared  to 
solving  the  Navier-Stokes  equations  over  the  entire  domain.  Equation  adap¬ 
tation  was  demonstrated  on  a  hybrid  sphere  mesh.  The  CPU  time  of  the 
viscous  tetrahedra  computation  was  67  %  greater  than  the  analogous  compu¬ 
tation  with  inviscid  tetrahedra.  This  is  a  significant  savings  in  computation 
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time.  A  further  advantage  to  inviscid  tetrahedra  is  the  memory  saved  by  not 
having  to  specify  extra  pointers  and  area  projections  of  the  tetrahedral  faces 
which  are  needed  for  viscous  calculations. 

6.  Use  of  spatial  adaptation  can  yield  accuracy  that  is  comparable  to  the  equivalent 
globally  fine  mesh .  Prismatic  spatial  adaptation  was  demonstrated  with  a 
three-dimensional  cavity  flow,  while  hybrid  adaptation  was  demonstrated  with 
sphere  flow.  Both  test  cases  illustrated  the  locally  adapted  and  globally  fine 
meshes  produced  similar  results. 

7.  Use  of  spatial  adaptation  provides  efficiency  by  yielding  substantial  CPU  time 
savings  relative  to  the  equivalent  globally  fine  mesh .  Prismatic  spatial  adap¬ 
tation  was  demonstrated  with  a  three-dimensional  cavity  flow,  while  hybrid 
adaptation  was  demonstrated  with  sphere  flow.  Both  test  cases  illustrated  the 
adapted  mesh  produced  the  same  solution  in  significantly  less  CPU  time  than 
the  globally  fine  mesh.  Memory  requirements  are  also  reduced  when  there  are 
fewer  nodes  or  cells. 

8.  The  developed  Navier- Stokes  scheme  is  robust .  However ,  it  could  be  pro¬ 
hibitively  expensive  for  periodic  time-accurate  computations.  The  scheme  con¬ 
verges  well  for  the  stretched  and  skewed  unstructured  grids  as  evidenced  by 
the  test  cases  on  adapted  and  hybrid  meshes.  However,  a  limitation  of  the 
scheme  lies  in  the  small  CFL  factor  necessary  for  stability  of  this  explicit  time 
marching  scheme. 

9.  The  Poisson  matrix  is  memory  intensive .  Unstructured  grids  created  around 
complex  geometries  frequently  contain  nodes  with  a  large  number  of  neigh¬ 
boring  points.  This  causes  a  substantial  increase  in  memory  requirements  for 
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the  storage  of  the  matrix  coefficients  of  the  Poisson  equation. 

7.4  Recommendations  for  Further  Research 

Based  on  the  present  work  on  simulating  viscous  incompressible  flows,  the  following 
recommendations  are  given  for  future  development. 

1.  Reduction  in  computation  time  would  allow  the  scheme  to  solve  a  wider  range 

of  steady  and  unsteady  problems  at  higher  Reynolds  numbers. 

(a)  With  explicit  time  marching,  the  time  step  is  restricted  by  the  CFL 
condition.  Implicit  schemes  are  not  restricted  by  the  CFL  condition  and 
larger  time  steps  can  be  employed.  A  semi-implicit  scheme  may  also 
remove  the  severe  CFL  limitation. 

(b)  A  comprehensive  comparison  of  several  matrix  solvers  (ICCG,  GMRES) 
may  produce  substantial  savings  in  computation  time. 

(c)  Application  of  a  multigrid  technique  for  the  Poisson  equation  could  be 
used  to  accelerate  the  convergence. 

(d)  Computing  requirements  may  be  reduced  by  performing  several  momen¬ 
tum  time  steps  before  correcting  the  velocity  back  to  a  divergence  free 
state  by  solving  the  pressure  Poisson  equation. 

(e)  Application  of  the  solver  on  parallel  computers  would  decrease  overall 
time  to  reach  a  solution. 

(f)  Local  time-stepping  would  allow  larger  time  steps  and  may  produce  faster 
convergence  for  steady  flow  computations. 
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2.  A  tradeoff  between  memory  and  CPU  time  can  be  examined  for  the  generation 
and  solution  of  the  Poisson  equation.  Currently  the  Poisson  matrix  entries  are 
calculated  initially  and  stored  in  arrays.  Since  the  finite  element  generation  of 
the  Poisson  matrix  is  relatively  simple,  calculating  the  matrix  coefficients  in 
the  iterative  solver  procedure  will  consume  more  CPU,  but  will  reduce  mem¬ 
ory  requirements  permitting  flow  simulations  about  larger  and  more  complex 
geometries. 

3.  Addition  of  a  turbulence  model  will  permit  the  application  of  the  method  to 
higher  Reynolds  numbers  flows. 

4.  The  non-zero  discrete  divergence  issue  needs  further  investigation. 
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Appendix  A 


Division  of  Prism  Cells  into 
Tetrahedra  Cells 


Each  prism  is  divided  into  three  tetrahedral  elements  using  the  method  of  [34].  To 
have  a  conservative  scheme,  the  tetrahedral  cells  created  from  the  prismatic  cells 
must  be  consistent  with  neighboring  tetrahedral  cells.  To  ensure  consistency,  the 
method  of  [34]  assigns  a  split  orientation  to  each  quadrilateral  face  of  the  prism,  as 
shown  in  Figure  A.l.  These  orientations  are  represented  by  arrows  on  the  upper 
triangular  face.  There  are  eight  possible  orientations  of  the  arrows,  but  two  of  those 


Split  Orientation 


Figure  A.l:  Prism  cell  division  showing  split  orientation  method. 
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orientations  do  not  result  in  division  of  the  prism.  These  orientations  are  when  the 
arrows  are  joined  nose  to  tail  in  a  clockwise  or  counter-clockwise  direction  around 
the  triangular  face.  The  division  of  prisms  must  be  completed  consistently  and  the 
two  non-permissible  orientations  must  be  avoided.  The  scheme  developed  by  [34] 
is  described  below: 

1.  For  each  node,  find  all  the  edges  containing  the  node.  Assign  an  orientation 
to  all  these  edges.  This  direction  points  to  the  current  node. 

2.  If  an  orientation  has  been  assigned,  take  no  action. 

This  scheme  will  guarantee  only  permissible  orientations  are  assigned.  For  the 
present  work,  the  first  layer  of  prisms  are  divided  using  this  split  orientation  method. 
Each  layer  of  prisms  repeats  the  pattern  of  division  developed  for  the  first  layer. 
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Appendix  B 


Data  Structures  for  Data  Arrays 

B.l  Primitive  Variable  Pointers 

Primitive  state  variables  and  pressure  values  at  the  grid  points  are  represented  by 
the  following  pointers. 


U (node,i),  i=l,3 

state  variables  corresponding  to  x,  y,  and  z 

momentum  equations 

U (node,i),  i=4,6 

change-in-time  of  the  state  variables  corre¬ 
sponding  to  x,  y,  and  z  momentum  equations 

P  {node) 

pressure 

DTNOD(  node) 

time  step 

VAR_N  (node) 

value  of  the  u  state  variable  at  the  previous 

time  step 

Arrays  for  the  pressure  correction  portion  of  the  solver  are  extensive.  The 
Poisson  matrix  is  stored  in  five  arrays,  one  array  holds  the  diagonal  terms,  a  second 
array  holds  the  terms  in  the  lower  diagonal  belonging  to  nodes  in  the  prism  region, 
and  a  third  array  holds  the  terms  in  the  upper  diagonal  belonging  to  nodes  in  the 
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prism  region.  Two  more  arrays  hold  the  terms  in  the  lower  and  upper  diagonal 
belonging  to  nodes  in  the  tetrahedral  and  interface  region.  Although  the  Gauss- 
Seidel  iterative  scheme  does  not  require  this  data  architecture,  it  would  permit  easy 
transition  to  a  more  advanced  matrix  solver.  Additional  arrays  are  necessary  to  hold 
the  neighboring  node  information  for  the  matrix.  The  divergence  of  the  velocity  and 
the  solution  vector  are  additional  arrays.  The  splitting  of  the  arrays  into  prism  and 
tetrahedral  regions  was  used  to  save  on  memory  storage.  The  number  of  neighboring 
nodes  for  a  node  in  the  prism  region  is  usually  less  than  20.  However,  due  to  the 
unstructuredness  of  the  tetrahedral  region,  the  number  of  neighboring  nodes  for  a 
node  in  the  tetrahedral  region  can  be  larger  than  80.  The  2-D  arrays  for  nodes 
in  the  prism  region  can  be  dimensioned  much  smaller  than  the  corresponding  2-D 
arrays  in  the  tetrahedral  region. 

DIAG (node)  diagonal  terms  of  the  Poisson  equation  matrix 

PHI  (node)  solution  vector  of  the  Poisson  equation 

DIVNODE  (node)  divergence  of  the  velocity  vector. 

INUM_U( node)  quantity  of  neighbor  nodes  of  primary  node. 

These  neighbor  nodes  are  greater  than  the  pri¬ 
mary  node 

INUMJL(node)  quantity  of  neighbor  nodes  of  primary  node. 

These  neighbor  nodes  are  lower  than  the  pri¬ 
mary  node 
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INDEX.U  ( i,  nodepr{sm ) 

INDEX-UT  (t,n«fetet) 

INDEX_L(  i,  nodeprism ) 

INDEX  ^LT(i,nodetet) 

AIC  C  G  _U  ( i,  nodeprism ) 

AICCG.UT  (i,nodetet) 

AIC  C  G  _L  ( i,  nodepriS7Tl ) 

AIC  C  G  _LT  ( i ,  nodetet ) 


node  neighbor  number  of  primary  node  in 
prism  region.  These  neighbor  nodes  are 
greater  than  the  primary  node.  The  index 
i  must  be  equal  to  or  less  than  the  number 
in  INUM_U( node) 

similar  to  INDEX_U(i, node),  except  pri¬ 
mary  node  is  in  tetrahedral  region, 
node  neighbor  number  of  primary  node  in 
prism  region.  These  neighbor  nodes  are 
lower  than  the  primary  node.  The  index  i 
must  be  equal  to  or  less  than  the  number 
in  INUM.L( node) 

similar  to  INDEX_L(i, node),  except  pri¬ 
mary  node  is  in  tetrahedral  region. 

Poisson  matrix  entry  of  that  neighbor  node 
stored  in  the  INDEX_U(i,  node)  location  for 
the  primary  node  in  prism  region, 
similar  to  AICCG_U(i; node),  except  pri¬ 
mary  node  is  in  tetrahedral  region. 

Poisson  matrix  entry  of  that  neighbor  node 
stored  in  the  INDEX_L(z; node)  location  for 
the  primary  node  in  prism  region, 
similar  to  AICCG_L(i,  node),  except  pri¬ 
mary  node  is  in  tetrahedral  region. 
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B.2  Node  Pointers 


Node  pointers  relate  properties  associated  with  each  node  in  the  mesh.  This  pointer 

is  oblivious  to  whether  the  node  lies  in  the  prismatic  or  tetrahedral  region. 

X  YZ(  0:node ,  i ) ,  i= 1 ,3  xyz  coordinates  of  each 

node.  i=l  refers  to 
the  x  coordinate,  i—2 
refers  to  the  y  coordi¬ 
nate,  and  i=3  refers  to 
the  z  coordinate 

VOLNOD(0:noc/e)  dual  cell  volume  of 

each  node 

ABSAREAJ3UALJ^OD(z,0:noe/e),  i=l,3  absolute  value  of  the 

xyz  area  face  projec¬ 
tions  for  the  dual  cell 
associated  with  each 
node.  This  metric  is 
used  only  in  the  time 
step  calculation 

B.3  Cell  Pointers 

Cell  pointers  relate  properties  associated  with  each  cell  in  the  mesh.  This  pointer  is 

oblivious  to  whether  the  cell  is  a  prism  or  tetrahedra. 

VOLCEL(0;ce//)  cell  volume  associated  with  each  cell 
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B.4  Edge  Pointers 

Edge  pointers  relate  properties  associated  with  each  edge  in  the  mesh.  These  point 


ers  are  oblivions  to  whether  the  edge 
AREAJEBGE(i,0:edge),  i=l,I 

IED  G  N  0  D  ( z, 0: edge) ,  i=l,2 
IUPD(0:ec/<7e,z),  i=l,3 

FC V(0:edge,i),  i=l,3 

GCV  (0:  edge, i),  i=l,3 

HCV(0:edge,i ),  i=l,3 


is  the  prismatic  or  tetrahedral  region, 
node-centered  dual  cell  face  xyz  area 

projections  associated  with  each  edge 
nodes  associated  with  each  edge 
working  array  used  in  the  coloring  the 
the  edges  and  nodes 
working  array  used  in  the  the  flux 
calculations 

working  array  used  in  the  the  flux 
calculations 

working  array  used  in  the  the  flux 
calculations 
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Appendix  C 


Data  Structures  for  Prismatic 

Grid 


The  prismatic  grid  consists  of  triangular  faces  that  cover  the  body  surface,  while 
quadrilateral  faces  extend  in  the  direction  normal  to  the  surface,  thus  constituting 
successive  layers  of  cells.  The  prismatic  grid  requires  a  set  of  pointers  to  define  a 
two  dimensional  triangular  mesh  combined  with  a  single  index  for  each  prismatic 
cell  that  belongs  to  the  same  stack  of  prisms.  Since  each  stack  of  prisms  is  formed 
from  a  triangular  face  on  the  body  surface,  a  prismatic  cell  may  be  referenced  by  the 
triangular  face  number  and  an  index  indicating  the  layer  of  which  the  cell  constitutes 
a  part. 


IBFACTOT 

IBNODTOT 

IBEDGTOT 

ICELTOT 

ILAYTOT 


number  of  faces  on  the  body 
number  of  nodes  on  the  body 
number  of  edges  on  the  body 
total  number  of  prism  cells 
number  of  prismatic  layers 
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NODTOT 


total  number  of  nodes  in  the  mesh 


IEDGTOT_PRISM  total  number  of  edges  in  the  mesh 

C.l  Face  Pointers 


The  triangulation  on  the  boundary  surface  that  constitutes  the  building  block  for 

the  prismatic  grid  is  represented  by  the  following  sets  of  pointers. 

IBFACNOD-TRI (ijace),  i=l,3  three  surface  nodes  associated 

with  the  surface  triangular 
face 

Using  the  semi-structure  of  the  grid,  the  nodes  that  constitute  a  face  on  layer 
i,  Figure  C.l,  are  given  by: 

N 1  =  NBl  +  (i-  1)*  IBNODTOT  (C.l) 

N2  =  NB2  +  (i  -  1)  *  IBNODTOT 
N3  =  NB3  +  (t  -  1) *  IBNODTOT 

C.2  Edge  Pointers 


The  prism  edges  can  be  numbered  locally  using  the  semi-structure  of  the  grid. 
IEDGTOT-VERT  is  the  number  of  vertical  edges  in  the  prismatic  region.  IEDG¬ 
TOT-VERT  is  equal  to  IBNODTOT*NL AYERS.  The  edges  that  constitute  a  cell 
on  layer  i,  Figure  C.l  are  given  by: 


El  =  EB1  +  IEDGTOT-VERT  +  (i  -1)*  IBEDGTOT  (C.2) 
E2  =  EB2  +  IEDGTOT-VERT  +  (i  -  1)  *  IBEDGTOT 
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N5 

E8 

N2 


NB1 

Figure  C.l:  Node  and  edge  numbering  system  for  prism  cells. 

E3  =  EB3  +  IEDGT0T-VERT  +  (i-  1)*IBEDGT0T 
EA  =  EBl  +  IEDGTOT-VERT  +  (i)*  IBEDGTOT 
E5  =  EB2  +  IEDGTOTJ/ERT  +  (i)  *  IBEDGTOT 
EG  =  EB3  +  IEDGTOTJ/ERT  +  (i)  *  IBEDGTOT 
El  -  N 1 
E8  =  N2 
E9  =  N3 

Note,  these  are  edge  numbers  based  on  the  local  edge  numbering.  The  global  edge 
numbers  obtained  by  hashing  the  edges  will  be  different  than  the  local  numbering 
based  on  the  prismatic  structure.  To  relate  the  two  numbering  systems,  the  pointer 
IEDG_GLOBAL(0;ed</e)  is  used:  given  a  local  edge  number,  the  global  edge  number 
is  produced.  This  is  not  a  necessary  pointer,  but  it  is  used  in  a  tradeoff  for  saving 
computation  time  at  the  expense  of  memory  storage. 
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N2 


Figure  C.2:  Boundary  edge-based  node  and  prismatic  cell  numbering  system. 

The  pointers  associated  with  the  edges  on  the  boundary  surface  are  illus¬ 
trated  in  Figure  C.2.  The  nodes  of  each  edge  on  the  boundary  surface  are  repre¬ 
sented  by  IB  ED  G  N  0  D  ( i,  hedge) ,  i=l,2.  The  neighboring  faces  (and  cells)  of  the  edge 
are  given  by  IBEDG_CELLS( t, hedge) ,  i=l,2. 


IBEDGNOD(l,IE)  =  N 1 
IBEDGN0D(2,  IE)  =  N2 
IBEDG-CELLS(l,  IE)  =  A 
IBEDG.CELLS{2,  IE)  =  B 


(C.3) 
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Appendix  D 

Data  Structures  for  Tetrahedral 

Grid 


The  cell  and  face  counters  for  the  tetrahedral  portion  of  the  grid  are  listed  below. 

The  nodes  and  edges  of  the  tetrahedral  cell  are  numbered  as  shown  in  Figure  D.l. 
ITCELTOT  number  of  tetrahedra  cells 

IEDGTOT  total  number  of  edges  in  the  mesh 

ITCELTOT JDIV  number  of  tetrahedra  cells,  including  the  di¬ 

vided  prisms 

IFSYMTOT_TET  number  of  faces  on  the  symmetry  boundary 

IFFFDTOT.TET  number  of  faces  on  the  far  field  boundary 

IFACTET-TOT  number  of  tetrahedral  faces.  Nonzero  if  vis¬ 

cous  tetrahedra  are  used. 
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N4 


Figure  D.l:  Node  and  edge  numbering  system  for  tetrahedral  cells. 


D.l  Cell  Arrays 


The  tetrahedral  portion  of  the  hybrid  mesh  requires  pointers  relating  the  nodes  with 
their  corresponding  cell.  The  volume  of  each  tetrahedra  is  also  stored  in  its  array. 
Not  only  are  these  nodes  and  volumes  stored  for  the  tetrahedral  portion  of  the  mesh, 
but  they  are  also  stored  for  the  tetrahedral  cells  created  by  the  division  of  the  prisms 

for  the  finite  element  discretization  of  the  Poisson  equation. 

ICELNOD(  Z,ceZZ),  i=l,4  four  nodes  associated  with  the  tetrahedra 

VOLTET (cell)  volume  of  the  tetrahedra 

D.2  Face  Arrays 

Cell-to-face  and  edge-to-face  pointers  are  required  for  the  case  of  viscous  tetrahedra. 
If  inviscid  tetrahedra  are  used  IFACTETJTOT  is  equal  to  zero,  although  there  are 
tetrahedral  faces  present. 
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ITFAC-CEL (0:face)  primary  cell  associated  with  a  face 

ITFAC_NOD( ijace),  i=l,3  three  nodes  of  the  face 
ITFAC_EDG(z,/ace),  i=l,6  the  first  three  indices  correspond  to 

the  three  edges,  not  on  the  face,  of  the 
primary  cell.  The  second  three  indices 
correspond  to  the  three  edges,  not  on 
the  face,  of  the  secondary  cell 

AREA_TET_FAC(i,/ace),  xyz  area  projections  of  the  face 

i=l,3 

D.3  Interface  Arrays 

A  few  arrays  are  necessary  to  define  the  pointers  needed  in  the  interface  region 
between  the  prisms  and  the  tetrahedra.  These  pointers  are  required  only  if  viscous 
tetrahedra  are  used. 

IBFAC_TCEL(/ace)  tetrahedral  cell  number  associated 

with  the  prismatic  boundary  face 

ITFAC  _EDG(z,/ace),  i=l,3  indices  which  correspond  to  the  three 

edges,  not  on  the  face,  of  the  tetrahe¬ 
dra  cell  on  the  interface. 
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Appendix  E 


Study  of  Unsteady  2-D  Flows 


Flow  around  circular  cylinders  is  of  significant  engineering  interest.  The  corre¬ 
sponding  flow  fields  are  quite  complex,  involving  such  flow  phenomena  as  boundary 
layer  separation,  vortex  formation  and  shedding,  as  well  as  vortex/vortex  interac¬ 
tion.  Study  of  two-dimensional  oscillatory  flows  is  of  importance  to  the  case  of 
wave-induced  forces  on  cylindrical  structures.  It  is  a  first  step  to  understanding  the 
complex  three-dimensional  wave- structure  interaction. 

The  present  study  applies  an  existing  two-dimensional  incompressible  Navier- 
Stokes  solver  to  examine  oscillating  flows  around  cylinders.  The  flow  will  be  defined 
and  results  of  the  two-dimensional  study  are  presented. 

E.l  2-D  Navier-Stokes  Solver 

In  the  present  work,  a  finite- volume/finite-element  numerical  scheme  has  been  used 
to  solve  the  unsteady  Navier-Stokes  equations  of  incompressible  flow  in  two  dimen¬ 
sions.  Quadrilateral  elements  are  used.  The  momentum  equations,  combined  with 
a  Simplified  Marker  and  Cell  (SMAC)  [4,  5]  type  of  pressure  correction  equation, 
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are  solved  employing  a  non-staggered  grid  where  all  of  the  dependent  variables  are 
defined  at  the  cell  corners.  The  momentum  equations  are  solved  using  a  finite- 
volume  numerical  scheme,  while  the  pressure  correction  equation  is  solved  using 
a  finite-element  numerical  scheme.  The  solution  is  advanced  in  time  with  an  ex¬ 
plicit/implicit  marching  scheme.  The  numerical  code  has  been  validated  in  a  previ¬ 
ous  study  for  a  flat  plate,  steady  uniform  flow  around  a  circular  cylinder  at  Reynolds 
numbers  of  40,  and  uniform  flow  around  a  circular  cylinder  at  Reynolds  numbers  of 
16200  and  105  [64,  91].  Further  details  on  the  2-D  solver  are  included  in  Appendix 
G. 


E.2  Forces  on  the  Cylinder  in  Oscillating  Flow 


Experimental  and  numerical  results  of  a  cylinder  in  oscillating  flow  show  the  in-line 
and  lift  coefficients  have  a  periodic  nature,  corresponding  to  the  periodicity  of  the 
freestream  oscillatory  flow,  which  has  the  form  of  U0 0  =  U0cosut.  Time  history  of 
the  in-line  and  lift  coefficients  are  usually  phase  averaged  so  that  the  in-line  and  lift 
force  during  a  single  averaged  period  can  be  examined.  In  the  present  work,  the 
numerical  simulations  were  run  for  a  minimum  of  ten  periods,  while  the  results  for 
periods  six  through  ten  were  phase  averaged.  The  time  history  results  show  steady 
periodic  results  after  the  fifth  period. 

In  the  present  work,  the  in-line  and  lift  force  coefficients,  Cx  and  Cl  respec¬ 
tively,  are  defined  as  follows 


Cx 


XF 

\PUID 


Cl  = 


YF 

\PVID 


(E.l) 

(E.2) 


where  XF  is  the  in-line  force,  YF  is  the  lift  force,  p  is  the  density,  f/^  is  the  freestream 
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speed,  and  D  is  the  cylinder  diameter.  The  periodic  in-line  force  is  dominated  by 
the  freestream  oscillating  flow,  whereas  the  periodic  lift  force  is  sensitive  to  the 
movements  of  the  vortices  as  they  separate  from  the  cylinder  and  are  swept  around 
the  cylinder  when  the  flow  reverses.  As  a  rough  rule,  the  fundamental  lift  frequency 
is  (m+1)  times  the  oscillation  frequency,  where  m  is  the  number  of  vortices  shed  per 
half  cycle  [142].  Therefore,  the  lift  force  becomes  a  strong  indicator  of  the  particular 
flow  regime  the  cylinder  is  experiencing. 

Numerical  results  from  the  two-dimensional  code  will  be  compared  against 
experimental  results.  Although  the  experiments  were  designed  to  be  two-dimensional, 
in  actuality,  the  shed  vortex  sheet  is  irregular  and  has  three-dimensional  coherent 
structures  superimposed  on  the  two-dimensional  coherent  structures  [53].  The  two- 
dimensional  numerical  simulation  cannot  resolve  any  three-dimensional  behavior, 
however,  the  major  characteristics  of  the  in-line  and  lift  force  curves  should  compare 
reasonably  well.  Minor  deviations  of  the  numerical  results  from  the  experimental 
results  is  expected. 

In  the  following  sections,  the  unsteady  flow  simulations  are  presented,  and 
comparisons  with  experimental  data  and  other  numerical  results  are  performed. 
Both  oscillating  flow  simulations  were  calculated  on  a  CRAY  C90. 

E.3  Instability  and  Numerical  Simulation 

Unsteadiness  in  the  form  of  the  Karman  vortex  street  occurs  during  any  physical 
experiment  of  flow  around  a  circular  cylinder.  Conceptually,  the  origin  of  the  de¬ 
struction  of  the  symmetric  pattern  can  be  explained  by  the  presence  of  multiple 
sources  of  perturbation  in  the  physical  model.  If  the  Reynolds  number  is  less  than 
40,  those  effects  are  dissipated  by  the  viscous  stresses.  However,  at  higher  Reynolds 
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numbers,  they  are  not  damped  and  the  flow  is  unsteady. 

In  the  case  of  the  numerical  simulation  of  unsteady  flow  past  a  circular  cylin¬ 
der  all  these  destabilizing  effects  are  absent.  Since  the  geometry  of  the  flow  and  the 
initial  and  the  boundary  conditions  are  symmetric,  the  simulation  of  the  Navier- 
Stokes  equations  leads  to  a  symmetric  solution  even  for  values  of  Reynolds  number 
greater  than  40.  The  truncation  and  round-off  errors,  as  well  as  those  due  to  the 
numerical  approximation  schemes,  are  perturbation  factors  which  could  eventually 
generate  vortex  shedding  [20,  21,  58,  59].  It  is  well  known  that  the  vortex  shedding 
can  be  generated  in  the  case  of  high  Reynolds  number  flow  if  there  are  no  or  lit¬ 
tle  explicit  /implicit  dissipation  terms.  However,  in  such  cases,  the  solution  will  be 
numerically  unstable. 

The  most  reasonable  way  to  generate  vortex  shedding  should  be  to  introduce 
into  the  numerical  model  the  same  perturbations  that  occur  during  a  physical  ex¬ 
periment,  if  these  perturbations  were  accurately  prescribed.  Unfortunately,  it  is  not 
possible  to  know  the  characteristics  of  these  perturbations  in  detail.  In  numerical 
simulations,  an  asymmetry  is  usually  introduced  into  the  initial  solution  in  order  to 
destroy  the  symmetric  flow  pattern.  Different  initial  movements  of  the  cylinder  have 
been  employed  in  the  literature  in  order  to  introduce  the  asymmetry  [20,  59,  89]. 
One  such  method  is  to  rotate  the  cylinder  for  a  short  time.  It  has  been  found  that 
the  final  periodic  flow  pattern  is  the  same  for  different  types  of  initial  perturbation 
of  numerical  simulations  [20].  This  fact  indicates  that  the  periodic  character  of  the 
flow  appearing  beyond  a  critical  Reynolds  number  is  an  intrinsic  property  of  the 
Navier- Stokes  equations  and  does  not  depend  on  the  nature  of  the  initial  perturba¬ 
tions.  These  conclusions  are  based  on  comparisons  made  by  varying  the  magnitude 
and  timing  of  the  perturbations.  It  was  not  reported  if  the  sign  of  the  perturbation 
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Figure  E.l:  Numerical  perturbation  to  induce  asymmetry. 


was  also  varied. 

In  the  present  work,  the  perturbation  introduced  corresponds  to  a  clockwise 
rotation  of  the  cylinder  followed  by  a  resting  period,  and  then  a  counterclockwise 
rotation,  shown  in  Fig  E.l.  Moving  the  surface  of  the  cylinder  causes  shear  stress 
in  the  flow  field  due  to  the  viscosity  of  the  fluid.  Using  this  kind  of  perturbation, 
vorticity  generated  through  shear  stress  in  the  region  around  the  cylinder  is  shed 
downstream.  During  the  perturbation,  one  side  of  the  surface  of  the  rotating  cylinder 
is  moving  in  the  direction  of  the  freestream  flow  and  a  lesser  amount  of  vorticity 
is  produced  compared  to  that  generated  on  the  other  side,  which  is  moving  against 
the  flow.  This  uneven  production  of  vorticity  destroys  the  symmetric  pattern  of  the 
flow  and  causes  unsteadiness  in  the  flow  field.  The  cylinder  is  slightly  rotated  only 
during  this  initial  perturbation.  During  the  rest  of  the  simulation  the  cylinder  is 
stationary. 
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E.4  Simulation  of  Oscillating  Flow  Around  Cylinders 


Oscillating  flows  around  cylindrical  bodies  are  quite  complex  and  significantly  dif¬ 
ferent  flow  patterns  can  be  created.  These  distinctive  flow  patterns  are  known  to 
be  a  function  of  the  Reynolds  Number  (Re)  and  Keulegan-Carpenter  Number  (KC) 
[118],  where  UQ  is  the  maximum  fluid  velocity,  T  is  the  period  of  the  oscillating  flow, 
Poo  is  the  freestream  density,  poo  is  the  freestream  viscosity  and  D  is  the  cylinder 
diameter. 


Re  =  PooU°D 

Moo 

(E.3) 

KC  =  V°T 

D 

(E.4) 

Often  the  frequency  parameter  is  used  to  categorize  oscillating  flows,  /?,  which  is 
the  product  of  KC  and  Re.  This  study  of  oscillating  flows  is  limited  to  two  distinct 
flow  regimes,  the  transverse  regime  and  the  double  pair  regime. 

E.4.1  Transverse  Flow  Regime,  KC  =  8,  Re  =  1568 

The  flow  regime,  when  the  KC  number  is  between  7  and  12,  is  called  a  transverse 
street  [142],  and  the  wake  is  comprised  of  the  shedding  of  one  large  vortex  during 
each  half  cycle.  The  transverse  street  is  composed  of  a  street  or  jet  of  vortices 
moving  away  from  the  cylinder,  roughly  perpendicular  to  the  oscillation  direction. 

The  lift  force  in  the  transverse  flow  regime  has  a  distinctive  shape,  which  is 
qualitatively  shown  in  Fig  E.2.  Since  there  is  one  vortex  shed  per  half  cycle  (m=l), 
the  fundamental  lift  frequency  is  twice  that  of  the  oscillation  frequency,  and  can  be 
identified  by  the  four  extrema  in  one  complete  period.  Each  negative  peak  in  the 
lift  force  is  caused  by  the  growth  and  shedding  of  a  large  vortex  in  each  half  cycle. 
Positive  peaks  are  induced  by  the  return  of  the  most  recently  shed  vortex  towards 
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Lift  Force 


One  Period 

Figure  E.2:  General  shape  of  lift  force  coefficient  for  transverse  mode. 

the  cylinder  just  after  flow  reversal  [142]. 

Numerical  parameters  used  for  this  case  of  KC  =  8,  and  Re  =  1568  are  shown 
in  Table  E.l.  The  far  field  boundary  is  located  at  approximately  33  diameters.  Per¬ 
turbation  parameters  used  to  induce  asymmetry  into  the  flow  are  shown  in  Table 
E.2.  Time  history  of  the  in-line  and  lift  force  coefficients  are  shown  in  Figures  E.3 
and  E.4.  Phase  averaged  in-line  and  lift  force  coefficients  are  compared  with  numer¬ 
ical  results  from  the  literature  [59].  These  comparisons  are  shown  in  Figures  E.5 
and  E.6.  The  phase  averaging  procedure  is  described  in  Appendix  F.  The  numeri¬ 
cal  results  from  the  literature  are  results  of  a  stream  function-  vorticity  calculation 
which  used  80  nodes  in  the  radial  direction  and  128  nodes  in  the  circumferential 
direction.  The  results  from  the  literature  were  for  one  instantaneous  period  and 
were  not  phase  averaged.  The  two  results  are  quite  similar  and  the  lift  force  shows 
the  flow  is  representative  of  the  transverse  mode. 
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Table  E.l:  Numerical  parameters  used  for  oscillating  cylinder  flow  at  KC=8  and  Re 
=  1568. 


Number 

of 

Radial 

Nodes 

Number 

of 

Circumferential 

Nodes 

Radial 

Stretching 

Factor 

Near  to 

Wall 

Step 

Size 

®  4u 

&4v 

&4w 

&4p 

96 

144 

1.08 

0.001 

1.0e-3 

5.0e-3 

Table  E.2:  Numerical  perturbation  parameters  used  for  oscillating  cylinder  flow  at 
KC=8  and  Re  =  1568. 
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Period 


Figure  E.3:  Numerical  in-line  force  coefficient  versus  time  for  transverse  mode, 
KC  =  8  and  Re  =  1568. 

- Freestream  velocity 

- In-line  force  coefficient 


Figure  E.4:  Numerical  lift  force  coefficient  versus  time  for  transverse  mode, 
KC  =  8  and  Re  =  1568. 

- Freestream  velocity 

- Lift  coefficient 
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Figure  E.5:  In-line  force  coefficient  comparison  for  transverse  mode,  KC  =  8 
and  Re  =  1568. 

o  Instantaneous  numerical  data  by  Justesen  [59] 

- Freestream  velocity 

- Phase  averaged  in-line  force  coefficient 


Figure  E.6:  Lift  force  coefficient  comparison  for  transverse  mode,  KC  =  8 
and  Re  =  1568. 

o  Instantaneous  numerical  data  by  Justesen  [59] 

- Freestream  velocity 

- Phase  averaged  lift  coefficient 
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One  Period 

Figure  E.7:  General  Shape  of  Lift  Force  Coefficient  for  Double  Pair  Mode 

E.4.2  Double  Pair  Flow  Regime,  KC  =  15,  Re  =  19050 

The  flow  regime,  when  the  KC  number  is  between  15  and  24,  is  called  a  double  pair 
[142],  and  the  wake  is  comprised  of  the  shedding  of  two  large  vortices  in  each  half 
cycle,  and  the  formation  of  two  pairs  in  each  cycle.  A  series  of  vortex  pairs  convect 
away  from  the  cylinder  in  opposite  directions  and  from  opposing  diagonal  points  on 
the  cylinder. 

The  lift  force  in  the  double  pair  flow  regime  has  a  distinctive  shape,  which 
is  qualitatively  shown  in  Fig  E.7.  Since  there  are  two  vortices  shed  per  half  cycle 
(m=2),  then  the  fundamental  lift  frequency  is  three  times  that  of  the  oscillation 
frequency,  and  can  be  identified  by  the  six  extrema  in  one  cycle. 

Numerical  parameters  used  for  the  case  of  KC  =  15  and  Re  =  19050  are 
shown  in  Table  E.3.  The  far  field  boundary  is  located  at  approximately  30  diameters. 
Perturbation  parameters  are  listed  in  Table  E.4.  Time  history  of  the  in-line  and  lift 
force  coefficients  are  shown  in  Figures  E.8  and  E.9. 

Obasaju,  et  al,[93]  have  found  mirror  image  modes  for  certain  flow  regimes. 
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Table  E.3:  Numerical  parameters  used  for  oscillating  cylinder  flow  at  KC=15  and 
Re  =  19050. 


Number  of 
Radial 
Nodes 

Number  of 
Circumferential 
Nodes 

Radial 

Stretching 

Factor 
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to  Wall 
Step  Size 

u 

®4v 

&4w 

®4  p 

in 

144 

i.i 

0.0002 

5.0e-3 

1.0e-3 

Table  E.4:  Numerical  perturbation  parameters  used  for  oscillating  cylinder  flow  at 
KC=15  and  Re  =  19050. 
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These  modes  and  their  corresponding  mirror  image  modes  both  generate  lift  of 
roughly  the  same  magnitude  but  of  opposite  signs.  Experimentally,  the  flow  ran* 
domly  switches  from  one  mode  to  the  other.  When  the  experimental  results  are 
phase  averaged,  the  lift  coefficient  curve  averages  to  zero.  To  obtain  useful  infor¬ 
mation,  Obasaju  separates  the  results  into  the  regular  mode  results  and  the  mirror 
image  mode  results.  Numerically,  the  random  switching  of  flow  modes  from  one  to 
the  other  does  not  occur.  However,  we  found  the  mirror  image  mode  can  be  pro¬ 
duced  by  reversing  the  numerical  perturbation  used  to  induce  asymmetry.  Results 
from  our  original  perturbation  technique,  shown  in  Figure  E.l,  corresponded  to 
Obasaju’s  mirror  image  mode.  Obasaju’s  original  mode  resulted  from  the  opposite 
numerical  perturbation,  a  counterclockwise  rotation  being  initiated  first,  followed  by 
a  resting  period,  and  then  a  clockwise  rotation  concluding  the  perturbation.  Phase 
averaged  in-line  and  lift  force  coefficients  are  compared  with  Obasaju’s  experimental 
results.  These  comparisons  are  shown  in  Figures  E.10  and  E.ll.  The  lift  force 
coefficient  is  the  mirror  image  mode,  so  we  have  plotted  our  results  against  the  neg¬ 
ative  of  Obasaju’s  result.  The  lift  coefficient  of  our  normal  mode  is  shown  in  Figure 
E.12  and  is  compared  against  the  actual  data  from  Obasaju.  The  numerical  results 
are  quite  similar  to  the  experimental  results  and  are  indicative  of  the  double  pair 
regime. 
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Figure  E.8:  Numerical  in-line  force  coefficient  versus  time  for  double  pair  mode,  KC 
=  15  and  Re  =  19050. 

- Freestream  velocity 

- In-line  force  coefficient 


Figure  E.9:  Numerical  lift  force  coefficient  versus  time  for  double  pair  mode,  KC  = 
15  and  Re  =  19050. 

- Freestream  velocity 

- Lift  coefficient 
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Figure  E.10:  In-line  force  coefficient  comparison  for  double  pair  mode,  KC  =  15 
and  Re  =  19050. 

o  Experimental  phase  averaged  data  by  Obasaju  [93] 

- Freestream  velocity 

- Phase  averaged  in-line  force  coefficient 
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Figure  E.ll:  Lift  force  coefficient  comparison  for  mirror  image  double  pair  mode, 
KC  =  15  and  Re  =  19050. 

o  Negative  of  experimental  phase  averaged  data  by  Obasaju  [93] 

- Freestream  velocity 

- Phase  averaged  lift  coefficient 


Figure  E.12:  Lift  force  coefficient  comparison  for  double  pair  mode,  KC  =  15 
and  Re  =  19050. 

o  Experimental  phase  averaged  data  by  Obasaju  [93] 

- Freestream  velocity 

- Phase  averaged  lift  coefficient 
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Appendix  F 


Phase  Averaging  Procedure  for 
Oscillating  Flow  Solutions 


The  time  varying  velocity  and  force  coefficient  solutions  for  oscillating  flows  around 
cylinders  are  quite  complex.  However,  the  solutions,  as  well  as  the  freestream  flow, 
may  have  a  periodic  nature.  After  several  periods,  the  velocities  and  force  coeffi¬ 
cients  often  reach  a  “steady  periodic”  solution.  Once  the  velocities  or  force  coeffi¬ 
cients  reach  this  steady  periodic  point  in  time,  the  results  can  be  phase  averaged . 
Phase  averaging  takes  the  solution  over  several  periods  and  averages  the  data  at  cor¬ 
responding  phases  of  the  period.  Once  the  time  varying  solution  is  phase  averaged, 
the  solution  can  be  displayed  for  a  single  averaged  period. 

F.l  Sampling  the  Time  Varying  Data 

The  phase  averaging  procedure  begins  with  sampling  of  the  data.  The  minimum 
time  at  which  the  phase  averaging  will  begin,  tmin,  and  the  maximum  time  at  which 
to  end  the  phase  averaging,  tmax ,  must  be  specified.  Examining  the  in-line  force 
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Figure  F.l:  Sampling  of  time  varying  data  within  specified  time  intervals. 

coefficient  at  KC=8,  Figure  E.3,  shows  this  solution  is  steady  and  periodic  after  the 
fifth  period.  Therefore,  sampling  from  the  sixth  to  the  twelfth  period  is  appropriate 
in  this  case.  Since  each  period  is  8  time  units,  tmin  =  48.0  and  tmax  =  96.0.  The 
time  increment  or  interval  for  the  integration,  tint ,  must  also  be  specified.  This 
interval  for  the  sampling  is  often  determined  by  numerical  experimentation,  but  a 
value  of  0.1  is  usually  a  good  choice. 

The  sampling  procedure  begins  at  the  value  of  time  =  tmin.  The  values  of 
all  data  points  which  fall  between  ti  =  tmin  and  t2  —  tmin  +  tint  are  averaged.  The 
value  of  time  in  this  interval  is  also  averaged.  The  multiple  data  points  within  this 
sampling  interval  are  replaced  by  the  averaged  single  data  point.  The  next  interval 
is  defined,  from  =  tmin  +  tint  and  t2  =  tmin  +  2*tint,  and  the  data  sampled. 
This  procedure  is  repeated  for  successive  sampling  intervals  until  t2  =  tmax.  This 
sampling  procedure  is  illustrated  in  Figure  F.l. 

The  averaging  procedure  in  the  sampling  process  is  actually  a  weighted  av- 
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Figure  F.2:  Weighted  averaging  of  sampled  data. 


eraging.  Each  data  point  is  weighted  with  its  own  associated  percent  increment  in 
time.  For  an  interval  i,  which  contains  the  data  points  a,b,c,  and  d,  the  sampled 
velocity  for  the  interval  is  shown  in  Figure  F.2.  The  final  sampled  velocity  for  the 
ith  interval  is  shown  in  Equation  F.l. 


y  /  -\  y  ^la  |  "h  Vb  tab  Vb  ~h  ^ 6c  .  Vc  "H  Vd  t cd  ,  y  t d2 

sample  a  tint  2  tint  2  tint  2  tint  dtint 


(F.l) 


F.2  Phase  Averaging  the  Sampled  Data 

Once  the  data  is  sampled,  the  phase  averaging  process  can  begin.  Since  the  sam¬ 
pling  time  interval,  tint,  as  well  as  the  period  is  known,  the  number  of  sampling 
time  intervals  within  one  period  of  time  can  be  determined.  The  sampled  data 
is  processed  by  averaging  the  entries  in  the  first  time  interval  for  all  the  sampled 
periods.  For  example,  if  there  are  5  periods,  with  10  sampling  intervals  in  each 
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period,  there  are  a  total  of  50  sampled  intervals.  The  data  at  intervals  1,  11,  21, 
31,  and  41  are  averaged  to  obtain  the  phase  averaged  data  in  the  first  interval  of 
the  phase  averaged  period.  This  averaging  is  continued  over  the  remaining  intervals 
in  the  period.  For  the  KC  —  8  example  mentioned  previously,  the  resulting  phase 
averaged  data  is  shown  in  Figure  E.5. 
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Appendix  G 


Description  of  2-D 
Navier-Stokes  Solver 


G.l  Governing  Equations  and  Pressure  Correction  For¬ 
mulation 

The  governing  equations  are  the  following  non-dimensionalized  laminar  Navier- 
Stokes  equations  of  incompressible  flow.  The  equations  were  non-dimensionalized 
using  the  freestream  velocity  and  the  cylinder  diameter. 


dV  dFi  5Gi_5Fv  dGy 
dt  +  dx  dy  dx  +  dy 


(G.l) 


The  state  vector  U,  the  convective  flux  vectors  Fi,  Gi,  the  viscous  flux  vectors 


Fv,Gv,  and  the  source  term  S  are  expressed  in  terms  of  the  primitive  variables, 


namely,  density  p,  and  x,y  velocities  and  pressure  P.  The  first  entry  in  the  vectors 


correspond  to  the  continuity  equation,  while  the  second  and  third  terms  correspond 
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to  the  x  and  y  momentum  equations,  respectively. 


An  explicit /implicit  marching  scheme  is  adopted  for  integration  in  time  of 
the  above  equations.  The  velocity  values  are  treated  explicitly,  while  the  pressure 
values  are  treated  implicitly  in  the  momentum  equations.  The  velocity  values  are 
marched  in  time  with  a  forward  Euler  scheme  [50].  The  continuity  equation  is 
formulated  implicitly  with  the  velocity  values  considered  at  time  level  ( n  +  1). 

The  momentum  equation  cannot  be  solved  directly  due  to  the  implicit  treat¬ 
ment  of  the  pressure  term.  An  auxiliary  velocity  vector  uf  is  introduced,  which 
satisfies  the  following  equation.  The  superscripts  denote  the  time  levels. 


dU/n+1  d*Y  dGf_  mr  dGy"  +1 

dt  dx  dy  dx  dy 


(G.2) 


In  the  auxiliary  velocity  vector  equation,  the  pressure  term  is  treated  explicitly 
and  u '  can  be  obtained  directly.  However,  the  solution  uf  does  not  satisfy  the 
continuity  equation.  Subtracting  the  auxiliary  velocity  vector  equation  (G.2)  from 
the  momentum  equation  (G.l),  it  is  obtained  : 


U<n+1)  -  t?  =  -  [v  At 


(G.3) 


Introducing  a  scalar  potential  such  that 


«<n+1)  -  u'  =  — V<£, 


(G.4) 
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the  following  equation  for  pressure  can  be  obtained: 

p(n+l)_p(n)  =  i_0  (G.5) 

Finally,  taking  the  divergence  of  each  side  of  equation  (G.4)  and  considering  the 
continuity  equation,  the  following  pressure  correction  Poisson  equation  is  obtained: 

A  <t>  =  -V-t?  (G.6) 

In  this  equation,  the  values  of  (f>  on  the  left  hand  side  are  treated  implicitly,  which 
requires  the  inversion  of  a  matrix.  Using  the  (j>  values  obtained  by  the  above  equa¬ 
tion,  we  can  correct  the  velocity  and  pressure  fields  using  equations  (G.4)  and  (G.5) 
as  follows: 


u<n+1)  =  u'  -V<f>  (G.7) 

P(n+1)  =  P(n)  +  (G.8) 

The  above  solution  procedure  follows  the  explicit /implicit  marching  scheme 
in  [71,  91].  The  overall  solution  procedure  corresponding  to  marching  by  one-time- 
step  is  summarized  as  follows: 

1.  Calculate  the  auxiliary  velocity  vector  uf  from  (G.2)  using  u(n)  and  p values. 

2.  Solve  the  pressure  correction  Poisson  equation  (G.6)  and  obtain  the  cj)  values. 

3.  Calculate  u (n+1)  and  p(n+1)  using  (G.7)  and  (G.8). 

4.  If  V  *  u(n+1)  <  €  where  e  is  the  tolerance  for  divergence,  advance  to  next  time 
step.  If  not,  consider  u(n+1)  as  uf  and  repeat  steps  2  and  3. 
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G.3  Finite-Element  Discretization  of  Pressure  Correc¬ 
tion  Equation 

The  pressure  correction  equations  are  discretized  using  the  Galerkin  finite-element 
approach.  The  scheme  is  compact  with  all  operations  being  restricted  to  within  each 
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grid-cell.  Bilinear  isoparametric  quadrilateral  elements  are  employed  [39]. 

For  quadrilaterals,  the  values  of  u,  v  and  <f>  are  defined  in  each  element  using 
the  following  finite-element  formulation  : 

4  4  4 

u  =  '£Ni-ui,  v  =  J2NiVU  <j>  =  Y,Ni&  (G.ll) 

1  =  1  1=1  1=1 

where  Ui,V{  and  &  are  nodal  values  of  u ,  v  and  <f)  and  N{  is  the  interpolating  shape 
function  associated  with  the  i-th  node. 

Integrating  equation  (  G.6)  over  the  each  element  domain  SI  using  the  Galerkin 
method,  the  following  equation  is  obtained  : 

Jj ^  Ni  {4>,xx  +  <t>,yy)  dxdy=  -  JJn  Ni  (u\x  +  v'y'j  dx  dy  (G.12) 

Applying  Gauss’s  theorem,  we  can  get  following  element  matrix  system  : 


(G.13) 


where 


Kx  = 

u 


D 


IL Ni^ dx  dy' K" = 


//.< 


dx 
dNi  dNi 


+ 


*3 

dNi  dN3 


q  \  dx  dx  dy  dy 


We  can  construct  the  following  global  matrix  system  by  assembling  the  element 
matrix  obtained  in  (  G.13) 


(G.14) 


where 


e 

f  =  2  [A'5 «;  +  K 5 »;] 
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The  values  for  on  the  left  hand  side  of  (  G.14)  are  treated  implicitly,  which 
requires  solution  of  a  system.  The  matrix  that  requires  inversion  is  a  symmetric 
linear  matrix.  The  system  is  solved  by  the  ICCG  ( Incomplete  Cholesky  Conjugate 
Gradient )  iterative  method  [94].  The  rate  of  convergence  of  the  iterative  process 
depends  on  the  numbering  of  nodes  [75].  In  order  to  accelerate  convergence,  the 
elements  of  the  matrix  are  renumbered  using  the  RCM  ( Reverse  Cuthitl-Mckee) 
method  [94]. 

A  typical  number  of  iterations  required  by  the  ICCG  method  for  the  Poisson 
equation  is  3,  while  5  iterations  are  typically  required  for  convergence  of  the  <f>  values. 
Convergence  criteria  for  the  pressure  corrections  is  V  <  10“2  'll/  L  where  U  and 
L  are  representative  velocity  and  length. 

G.3.1  Artificial  Dissipation 

Central  space  differencing  schemes  are  susceptible  to  oscillatory  modes  in  the  veloc¬ 
ity  field  of  high  Reynolds  number  flows.  Furthermore,  odd-even  decoupling  of  the 
solution  may  appear  in  the  pressure  field  for  this  non-staggered  type  of  mesh  that 
is  employed. 

In  the  present  work,  a  fourth  order  smoothing  term  is  added  explicitly  to  the 
momentum  equations  in  order  to  suppress  odd-even  decoupling  of  the  solution  [57]. 
Furthermore,  fourth  order  dissipation  is  added  to  the  pressure  correction  equation 
in  order  to  stabilize  the  solution  and  suppress  oscillations  in  the  pressure  field. 

The  smoothing  operator  is  cast  in  a  form  suitable  for  adaptive  unstructured 
grids.  All  operations  are  split  in  such  a  way  that  no  information  is  required  from 
outside  of  each  cell.  Each  grid  node  receives  contributions  from  each  one  of  its 
surrounding  cells.  The  operator  is  formed  in  two  steps.  The  second  order  difference 
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operator  is  formed  in  the  first  step.  The  second  order  distributions  to  cell-corners 
(j)  for  the  momentum  equations  are  as  follows  : 


D](u)  =  -  4uj 


(G.15) 


The  second  step  duplicates  the  first,  replacing  state  variables  by  second  order 
differences  from  the  first  step.  The  fourth  order  smoothing  distributions  are: 

-■£>;(«)  =  (l>,W)  - (G.16) 


\t  =  l  / 

The  fourth-order  difference  terms  of  u  and  v  are  multiplied  by  empirical  coefficients 
<J4(u)  and  ct4(v),  Large  values  of  a4  stabilizes  the  solution  but  destroys  the  accu¬ 
racy.  Therefore,  special  care  is  required  for  choosing  the  value  of  a4.  Numerical 
experiments  have  been  carried  out  to  determine  optimum  values  for  the  smooth¬ 
ing  coefficient.  The  determined  values  are  such  that  the  solution  accuracy  is  not 
affected,  while  the  odd-even  modes  are  suppressed  [61]. 

Similarly,  fourth  order  difference  terms  of  pressure  are  added  to  the  right 
hand  side  of  the  Poisson  equation  (  G.14)  as  follows  : 


D3>  =  f  +  d 


(G.17) 


where 

dT  =  [<74(p)Df  (pju))] 

Generally,  the  value  of  a4  (p)  is  different  from  a4  ( u )  and  a4  ( v ).  In  the  present  work, 
the  values  of  o4  (u)  and  a4  (u)  are  the  same,  while  the  value  of  a4  (p)  is  different. 


G.3.2  Time-Step  Calculation 

Using  central  space  and  forward  time  differencing,  the  stability  limitation  for  the 
model  1-D  convection  equation  ut  +  cuy  =  0  is  <  1  (CFL  limitation),  while  the 
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corresponding  stability  restriction  for  the  1-D  model  diffusion  equation  ut  =  vuyy 


;s  v*t  <  1 
-  2- 


In  the  present  scheme,  a  combination  of  the  two  limitations  is  employed. 
Specifically, 


At  =  umin{- 


Am 


A  l 


u\  + 


aim  M  +  S27 


}, 


(G.18) 


where  Am,  Al  are  the  cell  dimensions  in  the  m,Z  local  cell-directions,  u,v  are  the 
corresponding  velocity  components,  v  is  the  kinematic  viscosity  coefficient,  and  a  = 
\  is  the  diffusion  stability  limitation.  Lastly,  u  is  a  safety  factor,  and  equal  to  0.9. 


G.3.3  Boundary  Conditions 

Two  types  of  conditions  have  been  applied  for  the  cases  considered  in  the  present 
work:  wall  and  far  field.  The  boundary  conditions  are  applied  to  the  two  velocity 
components,  as  well  as  to  the  pressure  corrections.  At  the  wall,  the  u  and  v  com¬ 
ponents  of  velocity  are  set  to  zero.  The  value  of  d<f>/dn  is  also  set  to  zero.  At  the 
farfield ,  the  velocity  components  are  set  equal  to  the  free-stream  values,  while  the 
value  of  d(f)/dn  is  set  to  zero. 

An  added  complication  to  the  boundary  conditions  exist  because  the  pressure 
correction  equation  is  a  Poisson  equation,  and  derivative  boundary  conditions  are 
used.  In  this  case,  there  are  an  infinite  number  of  solutions  for  the  pressure  field 
differing  only  by  arbitrary  additive  constants.  In  these  circumstances  the  iterative 
solution  method  may  never  converge  unless  the  pressure  correction  <f  is  specified  at 
one  point  in  the  flow  [100].  In  the  present  work,  </>  at  a  single  upstream  point  on  the 
farfield  boundary  is  specified  to  be  zero. 
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